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ABSTRACT 


The objective of this thesis is to develop a technique 
and associated algorithms to extract the arrival time of modal 
energy, uSing a vertical array, from broadband signals. Modal 
energy arrival time is important to shallow water acoustic 
tomography because low angle rays, which contain the majority 
of acoustic energy, are often not resolvable. Tile 
compensation is included in the beamforming algorithm to 
provide a virtual vertical array. A broadband modal filtering 
technique is accomplished through weighting the frequency 
components of phase encoded tomographic signals by the 
spectrum of the mode shapes. A methodology of phase decoding 
after beamforming was adopted to minimize processing. Initial 
development and prototyping was done using a parabolic 
equation model. Further testing was accomplished on real data 
taken from the Barents Sea Polar Front Experiment, August 
1992. Results show consistency over a number of transmitted 
pulses. Mode energy travel time measurement is simplified due 
to the distinct arrival structure of beamformed signals. 
Based on these results, the modal beamforming algorithm should 


be a useful tool for acoustic tomograpy. 
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I. INTRODUCTION 


A. BACKGROUND 

Tomography signal processing requires the precise 
measurement of acoustic travel time, which is the integral 
along the raypath of inverse sound speed. Travel time of 
acoustic signals has been well determined to be a function of 
temperature, salinity, and pressure. Once time perturbations 
for multiple arrival paths are known, ocean fluctuations can 
be determined using mathematical inverse techniques [Ref. 1]. 
In addition to measuring time perturbations based on ray 
arrival times, it is possible to measure perturbations based 
on mode arrivals. There are two approaches to tomography 
using modes: for continuous wave signals, the basic 
measurement is the modal phase difference; for broad band 
signals, the measurement used is modal travel time [Ref. 2]. 
Since coded signals (maximal-length sequence phase encoded 
Signals) are intrinsic broadband signals, modal travel time 
measurement is germane to the work done in this paper. Where 
ray theory is not applicable, or not convenient, it is 
possible to extend tomographic techniques using both rays and 
modes [Ref. 3]. 

The goal of my thesis is to develop a technique for 


broadband modal beamforming for acoustic tomography. The 


scope entails broadband modal beamforming with a vertical 
array and decoding (m-Ssequence removal of phase encoded 
Signals). In the process of my research several software 
programs were developed. Most notable is the broadband modal 
beamformer, which is an adaptation of the continuous wave 
beamformer developed by Crocker [Ref. 4]. Specific goals of 
research include: 


® Develop a broad spectrum modal filtering technique given 
the normal mode’s eigen function and values. 


® Evaluate compatibility for phase sequence removal prior to 
and following beamforming. 


® Minimize processing time, where possible, and provide for 
a robust environment. 


@® Evaluate performance using both synthetic and in situ 
data. 


® Provide a virtual array based on, time varying array tilt, 
array geometry and modal group speed. 


® Incorporate a supporting cast of programs, including a 
modal decomposition program for standardizing input to the 
beamformer and to simplify and enhance future signal 
processing. 
B. THESIS SUMMARY 

Chapter II, Barents Sea Experiment, is a summary of the 
Barents Sea Polar Front Experiment (BSPFEX). This experiment 
is addressed because it provided real data to evaluate the 
processing algorithms. 


Chapter III, Acoustic Normal Mode Propagation Theory, is 


an introduction to normal modes. Modal acoustic propagation 


is addressed. The chapter is concluded with a section linking 
acoustic ray and normal mode theory. 

Chapter IV, Signal Processing Concepts, presents the 
concepts and theory behind broadband modal beamforming. 
Subjects include beam steering, mode filtering and phase 
encoded signals. A plane wave beam steering approach is taken 
to compensate for array tilt. Following tilt correction, my 
modal filtering methodology is presented. Finally, a brief 
introduction to m-sequence signal decoding is given. Phase 
preservation after beamforming is addressed, supporting the 
concept of decoding after beamforming. 

Chapter V, Barents Sea Array, is a description of the 
array used during the Barents Sea Polar Front Experiment 
(BSPFEX). Improvements made from previous vertical arrays are 
discussed. Statistics on hydrophone performance as well as 
the array modal beam patterns are addressed. 

Chapter VI, Evaluation of Barents Sea Data, is a 
presentation of results from the BSPFEX. Comparisons and 
conclusions are drawn concerning the array’s performance. 

Chapter VII, Summary and Future Work, presents a summation 
of conclusions and identifies areas requiring improvement. 

The appendix contains the most significant algorithms I 
developed in the course of research. Algorithms included are: 
the broad band modal beamformer, modal decomposition program 
and the supporting modal extraction program which formats the 


input files for use in the beamformer. 


II. BARENTS SEA EXPERIMENT 

Concepts developed in this thesis were applied to real 
data collected from the Barents Sea Polar Front Experiment 
(BSPFEX) . During the BSPFEX, I developed and tested my 
beamforming techniques and algorithms on synthetic data 
modeled using the Finite Element Parabolic Equation (PE) Model 
developed by Collins and Westwood [Ref. 5]. The synthetic 
data was modeled for the anticipated location of the array for 
the BSPFEX. Unfortunately, the array was deployed in a 
different location from the model limiting any comparisons 
drawn between the synthetic and BSPFEX data. 

The BSPFEX was conducted during August of 1992. The 
experiment comprised a joint effort between the Naval 
Postgraduate School (NPS), Woods Hole Oceanographic 
Institution (WHOI), and the Science Applications International 
Corporation (SAIC). The principal investigators for the 
experiment are Professors Robert Bourke, Ching-Sang Chiu, and 
James H. Miller from NPS, Dr. James F. Lynch and Dr. Al J. 
Plueddemann from WHOI, and Dr. Robin Muench from SAIC. The 
principal engineer for the vertical hydrophone array system 
was Mr. Keith Von der Heydt from WHOI. The objectives of the 
experiment as outlined by the Barents Sea Polar Front Group 


1992 are: 











alae Provide a detailed physical description of the polar 
masOrit. . 


Ze Enhance the understanding of dynamics of the front, 
including frontogenesis and its influence on regional 
oceanographic processes. 


or Assess the ability of acoustic tomography to define 
frontal and associated mesoscale features. 


4. Provide improved documentation of shallow water acoustic 
propagation in this region and the effect of the environment 
on acoustic ASW operations. 


[Ref. 6] 


This experiment is a milestone in acoustic tomography. It 
is the first time that a vertical array has been used in a 
shallow water environment to conduct acoustic tomography. A 
sixteen hydrophone telemetered vertical array, with ten meter 
spacing between hydrophones, was deployed during’ the 
experiment. The source and array mooring locations are listed 
in Table I. Table II describes the characteristics of the 
deployed broadband sources. 

Table I MOORING LOCATIONS FOR THE 224 AND 400 HZ SOURCES 


AND THE VERTICAL ARRAY. 


Mooring Latitude Longitude Depth (m) 


SE(a) VLA a oo ON eo 5 ooo 2158 

SE(b) VLA Wael LOG New 23°" 32.2960" 275: 

NE 224 Hz 1 eS 3 5 Nee 2 3 oa 2d SOS 1A! 

NW 400 Hz Va-eo2.,9Lo4 N 2° 44,1043" Ow 

SW 400 Hz 742 04 7337°N 22° 200. 46057 FE oe8 30% 

———E—E——E—E—E——— ee ————eE——eEanHRUEHHLHL_ESE_— SEE ee ee 
5 


© 0 @ © 


Table II BSPFEX SOURCE CHARACTERISTICS. 
aa eae ae A al 


Mooring NE SW NW 
Freq (Hz) 224 400 400 

SL (dB) 183 183 ilyers 
Duration (s) Li3.25 132.86 132.86 
BW (Hz) 16 100 100 

Q (cyc/dig) 14 4 4 
No. digits 63 Salk 5a 
Seq length 3.9305 Sra Lal 5 oat 
No. sequences 30 24 +2 24 +2 
Cycle (min) °2°57 75a ke 20. 5 aus, 25 


Overall, the experiment was a sterling success. The WHOI- 
NPS vertical array was deployed on 12 August 1992. The array 
performed flawlessly for fourteen hours when a failure in the 
amplifier chip forced retrieval and repair of the array. 
After repair, the array was redeployed without the adjustable 
gain control making signal quantization from analogue to 
digital less optimal. Following recovery and redeployment, 
the array operated continuously for three days recording 
signals from both the 224 Hz and 400 Hz sources. In ‘ane 
fifteen gigabytes of data were recorded. 

One major disappointment was the failure of the NW 400 Hz 
transceiver to transmit shortly after its deployment (see 
Figure 1). No attempt was made to recover and repair the 
broken transceiver “So aS to Hotesdi ctmmb the —tonoc maa noms 
deployment schedule, and in the hopes that the transceiver was 
Still able to receive acoustic data. Post experiment recovery 


of the broken transceiver revealed that it was still able to 





Figure 1 Deployment locations of the telemetered array 
and sources, from the BSPFEX. SW - 224 Hz transceiver; SE 
- 16 channel vertical array; NW/NE - 400 Hz transceivers 
iRet. 7]. 





receive and recorded three days of acoustic data. With the NW 
transceiver still partially operational, three paths for 
tomography are available; two across the polar front and one 
along the front. 

In addition to tomographic pulse receptions, other 
acoustic observations were conducted. The ship dropped eight 
SUS charges to measure reverberation and bottom properties. 
Eighteen additional SUS charges were dropped from a P-3 
aircraft along the southwestern and northern tracks, to be 
used in evaluating propagation and bottom loss. One hour of 


Continuous wave (CW) 224 Hz transmission was made for use in 


modal phase and matched field tomography, and to measure 
temporal coherence for comparison with pulsed transmissions. 

Other sensors used during the deployment include, over 
three hundred Conductivity-Temperature-Depth (CTD 
deployments, two of three current meters (one failure), and 
acoustic doppler current profiler (ADCP). Figure 2 denotes 
the fine structure of the polar front present, as measured 
from CTD deployments. Acoustic propagation was excellent 
during the experiment. The 224 Hz and 400 Hz sources were 
clearly audible at ranges over 50 kilometers, even without 
Signal processing. With such clear receptions and low noise 
levels, the task of processing the data later will be much 
simplified. The above summary was based on the preliminary 


cruise report for the BSPFEX [Ref. 7]. 


Temperature — Dense Section 2 (S->N) 


ASAE NORMAL IAIN 8 RAS AAA Or 


Nr 


km. 





Figure 2 Contour plot of the Barents Sea Polar Front 
Peer. 7). 


III. NORMAL MODE PROPAGATION THEORY 


A. ADIABATIC MODE THEORY 
Normal mode propagation theory is a physical model whose 


eigenfunctions, are solutions to the acoustic wave equation, 


0 





1 Op 
Wijoys = : (35) 
ee Ole 
in a waveguide [Ref. 8]. The objective of normal mode theory 


is to model the "local" pressure field as a linear combination 
of vertical standing waves. The term local refers to the 
range dependent sound speed structure and boundary conditions 
at a certain range. The following development applies fora 
slowly varying range dependent waveguide and is taken from 
Shang [Ref. 2]. 

The acoustic pressure field can be expressed in terms of 


local mode superposition, 


p(r,z) ya (3.2) 
n=1 


where W,(z,r) is the local solution or eigenfunction, and 
describes the distribution of modal energy as a function of 
depth. F,(r) 1s the range dependent portion of the pressure 
field, and the Vr term describes the geometric transmission 


losses caused by cylindrical spreading. 


10 


Subject to the boundary conditions and orthogonality, the 
local mode WV, must satisfy, 


3 1 OW, il . Pri . (3.3) 
Est) eal ea ibe = 0 


where k’(r,z) is the acoustic wave number, k,?(r) is the 
vertical wave number, or the eigen value, associated with the 
mode W,. p(z) is the depth dependent density. The range 


dependent portion of the pressure field must satisfy, 


aF 


dake i ' 
+ Ka(2 + ao \r ~ “2, [2A (7) PIL iG), En aa (3.4) 


r? 











Aw, and B,, are coefficients describing the first and second 


order coupling effects in the sound channel: 


E 1 OWE Zar ) 


Z a Oy (2Z,r) 3.6 
Ban = | S—gy Vn (212) —S az. (3.6) 


Tf mode coupling is weak enough, then coupling can be 


neglected, and we have the adiabatic solution, 


Jpatera' 
p(xr,z) = e*/*/any Yalelorvateley G7) 
m haa 


iL al 


B. NORMAL MODES AND ACOUSTIC TOMOGRAPHY 

The travel time of each normal mode is dependent on the 
ocean structure. It is possible to estimate mode travel time 
from acoustic models, or it can be directly measured using 
mode filtering techniques. The difference between measured 
and estimated arrival times are the time perturbations for 
each mode. Broadband modal acoustic tomography seeks to 
exploit these time perturbations by applying inverse 
techniques to infer information about the acoustic 
environment. 

For a broadband source, the pressure field can be 
described as a function of a continuous wave field modified by 


the source signal spectrum, 


oe 


p(r,zZ,t) = [s(o) pa (r,2,0) eF?*de, (me 


where s(o) is the signal spectrum of the source, and p,, is the 
continuous wave field. The total pressure is the sum of the 


pressures for every mode, 


Ditgaz,, Ce eae 2 oe (3.9) 


M 


where the modal pressure term is, 


jo Genes, f)) = [so BatFol Wal2) a 5tnte-00) ay, (3.10) 


se /K,,2 


eZ 





Taking the derivative of the phase with respect to frequency, 





to 2 (3.11) 
W 


and it follows that the group speed must be, 


Cr ® nae (3.12) 


m 

It has been shown that there iS a unique travel time 
associated with each mode for a given environment [Ref. 2]. 
Acoustic tomography can be accomplished uSing the time 
perturbations taken from the measured and estimated arrival 
times. The remaining task is to develop a technigue to 
measure the modal arrival times, which is the subject of the 


following chapter on signal processing. 


ae 


IV. SIGNAL PROCESSING CONCEPTS 


A. BROAD BAND MODAL BEAMFORMING 

Wave propagation modeling can be used to extract temporal 
and spatial characteristics for a propagating wave, and 
through inverse techniques to infer information about the 
acoustic medium. This is the essence of acoustic tomography. 
In recent years, ocean modeling for acoustic tomography has 
been dominated by acoustic ray theory. Ray theory is 
attractive because it is simple to model and requires minimal 
computer resources. Other modeling techniques simply were not 
practical until recent years. Ray theory is a high frequency 
approximation which fails to describe turning points, 
diffraction leakage and caustics. Further, many ray arrivals 
may not be resolvable, particularly in a shallow water 
environment. 

A vertical array provides an added dimension to 
tomography. Not only does it increase signal gain, it also 
allows mode and plane wave beamforming as well as enhancing 
resolvability of ray arrivals. The rest of this section will 


be devoted to the concepts required for modal beamforming. 
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1. Beam Steering 

Array tilt can be caused by a number of forces such as 
waves, ocean currents, and external tensioning. Any tilt in 
the array 1S undesirable, and must be corrected. A virtual 
vertical array can be established using a plane wave beam 
steering methodology. 

There are essentially two approaches to plane wave 
beam steering. The first applies phase shift delays to the 
Signal at the desired frequencies. The second uses true time 
delays. Phase delay beam steering is well suited for CW 
Signals. However, this method is much less efficient for 
broad band signals because it requires phase shifting for 
every frequency of the broad band signal. 

Time delay beam steering is more suitable for use with 
broadband signals. Time shift delays are determined from 
array geometry and sound speed. For a line array, the time 
shift delay for the n, element can be described as, 


_ ndsin(@,) 


= (4.1) 
E 


nN 


where @, is the tilt in the direction of the source, andd is 
the element spacing [Ref. 9]. For normal mode beamforming, c 
is the group speed, equation (3.16), for each mode. Equation 
(4.1) is elementary in nature, and can be applied to any array 


or signal. Applying shift delays to discrete time signals 


iS 


requires interpolation incurring additional processing [Ref. 
4). 
2. Mode Filtering 
Plane wave beamforming does not account for the non- 
uniform energy distribution in the wave guide as described by 
normal modes. There are two types of modal beamforming, 
continuous wave (CW) and broad band (BB). 

For CW modal beamforming, the task is simple. The 
time domain pressure field need only be weighted by the 
eigenfunction corresponding to each element depth. Phas 
method is not appropriate for broad band signals, because 
normal modes are a function of both depth and frequency. For 
BB beamforming, the frequency domain signal must be weighted 
by the corresponding frequency component of the eigenfunction. 
This requires that the mode shape be known at each frequency 
of the Discrete Fourier Transform (DFT) for the pressure field 
covering the desired bandwidth. 

For encoded signals with long duration sequence lengths, 
modal decomposition processing can be unsurmountable. The 
nominal modal decomposition time is thirteen minutes for each 
frequency. A 100 Hz bandwidth source with a sequence length 
of 5.11 seconds has 1022 frequency components. With current 
computer resources (SPARC 2), it would take over nine days to 
process the eigen functions for each frequency. Fortunately, 


the mode shapeS are continuous functions in depth and 
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frequency. Figure 3 shows how the mode amplitude, at a fixed 
depth, varies slowly with frequency. The mode shapes can be 
sampled as a function of frequency and then interpolated to 
achieve and desired frequency increment within the band. 
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Figure 3 Mode twenty as a function of frequency at 126 
meter fixed depth. 


The next step is to weight the frequency domain 
pressure field by the mode shape at depths corresponding to 
the array elements. The following equation is the essence of 


_ modal ea ae ering. 


a,(t) ge ct Pi Zee ee, f) ||, (4.3) 


diay, 


where a,” is proportional to the intensity of mode m, and p is 
the density in the water column. The only limitationmese 
equation (4.2) is that the Fourier Transforms are long enough 
to accommodate the length of the entire coded sequence. The 
pressure field frequency components are weighted by the 
frequency domain normal mode at the appropriate array element 
depths. The resultant is summed over depth, and an inverse 
Fourier Transform is applied to achieve the time domain 
beamformed signal. 

The array gain (AG) exhibited by broadband modal 
beamforming cannot be greater then conventional beamforming 


theory, 


AG = 10log(N), (4.3) 


[Ref. 10]. The theoretical array gain shown in equation (4.3) 
assumes the noise between elements is uncorrelated. Signal 


enhancement is strongly dependent on this correlation factor. 


B. MAXIMUM LENGTH PHASE ENCODED SIGNALS 

The objective of signal processing for broad spectrum 
acoustic tomography is the measurement of travel time for 
different raypaths and mode arrivals along a vertical cross 
section of the ocean. 


One method in which travel time can be measured is through 


the use of explosive and implosive devices. 
Unfortunately, explosive devices do not propagate 
repeatable acoustic signals. A better method that has 


been found employs the use of maximal-length sequences (m- 
sequences) or pseudo-random noise are well suited for this 
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application because of their deterministic nature, 
correlation properties, and simplicity [Ref. 11]. 


Beeequences are a sequence of binary bits having the property 
of a maximum possible period for an r-staged shift register. 
Shift registers are created from primitive polynomials. These 
binary shift registers are mapped to a series of ones and 
minus ones [Ref. 12]. The primary characteristic of the 
maximum length shift register is its autocorrelation 
properties. The self correlation of an m-sequence of length 
N produces a triangular function of short duration as shown at 


Figure 4. 





Figure 4 Autocorrelation function for a m-sequence of 
length N. 


An m-sequence coded signal is formed using a simple 
harmonic source and applying a phase shift from the m-sequence 


register. A phase encoded signal has the form, 


uo 


s(t) = Acos(2nf,t+M(t)o), (4.4) 


where A is the amplitude of the transmitted signal and M(t) is 
the value of the shift register for a specified digit 
duration. OQ is an integer number of cycles (at carrier 
frequency) per digit and is used to determine the digit 
duration. The phase modulation angle @ is chosen to optimize 


signal-to-noise performance: 


d@ = tan ?(/N). (4.5) 

The one major drawback to using coded signals is that the 
standard correlation requires N*’ multiplications. Decoding 
algorithms can over come this processing issue by using the 
Fast Hadamard Transform (FHT). The FHT is analogous to a FFT 
and reduces processing to NlogN multiplications. 

Decoding signals is a simple process. The first step is 
to filter the input signal using a band pass filter covering 
the frequencies containing the coded signal. After filtering 
the signal is demodulated to remove the carrier frequency. 
The demodulated signal is then correlated for different shift 
versions using the FHT to reduce processing. Figure 5 is a 
block diagram showing the Birdsall-Metzger (BM) decoding 


process used in [Ref. 11]. 
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Figure 5 Block diagram of m-sequence removal algorithm 
(Ref. 11] 
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C. SIGNAL PROCESSING SUMMARY 

The first step in modal beamforming is the correction of 
tilt to provide a virtual vertical array. This is 
accomplished using plane wave beamforming techniques. True 
time shift delays are used vice phase shift beamforming to 
reduce procesSing. After tilt correction, mode filtering is 
accomplished in the frequency domain. Modal beamforming 
constitutes filtering the frequency domain signals by the mode 
shapes, and then summing across the array elements. Post 
beamforming, the signals are decoded using the BM phase 
Gecoder algorithm. 

The phase decoding process can be accomplished prior to 
beamforming, but this would significantly increase signal 
processing as a multiple of the array elements and is not 
recommended. Further, the beamforming algorithm developed in 
this thesis does not Support this procedure. If at some point 
it becomes desirable to process signals in this manner, the 
beamformer program would have to be altered to except complex 
data and to perform mode filtering about zero frequency offset 


signals. 
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V. BARENTS SEA ARRAY 


A. ARRAY COMPONENTS AND DESIGN 

The WHOI-NPS array contains sixteen hydrophone elements at 
ten meter intervals. The telemetered vertical array was 
moored to the bottom uSing a 1000 lb. cast anchor. The 
distance from the anchor to the bottom hydrophone is 
approximately 3 meters. A forty-eight inch bouyant syntactic 
foam sphere suspends the array 166 meters above its mooring. 
One meter below the sphere are the dual Benthos Interrogators. 
The interrogators were used to acoustically determine the 
exact location of the tethered array and to measure array 
tilt. Two inclinometers, one located two meters beneath the 
interrogator and the second near the center of the array, 
provided an additional means of determining array tilt. 
Figure 6 is an illustration of the array design. 

Nominal hydrophone sensitivity was -160 dB re 1 v/yPa with 
a low frequency roll-off at 50 Hz. Data was preprocessed with 
an 8-pole Butterworth low pass filter at 500 Hz cutoff 
frequency. A sampling rate of 1600 Hz was selected to achieve 
twice Nyquist frequency of the highest frequency source 400 
Hz. A sampling frequency four times carrier frequency is 


required by the m-sequence removal algorithm to simplify 
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Figure 6 Diagram of the WHOI-NPS vertical array used during 
the Barents Sea Polar Front Experiment [Ref. 7]. 











Signal processing. Data from the 224 Hz source had to be 
resampled to 896 Hz sampling frequency for decoding. A single 
data bus was used to sample all sixteen array channels and 
created a 30 microsecond skew between elements. The 
beamformer corrects for the sampling skew based on the 


principles established for tilt compensation. 


B. BARENTS SEA POLAR FRONT EXPERIMENT - MODAL BEAMFORMER 
1. Equipment Performance 
The WHOI-NPS array was deployed in the Barents Sea on 


12 August 1992" The array performed well with one minor 
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mishap when the amplifier chip failed. The amplifier chip 
failure effected adjustable gain control and degraded analogue 
to digital quantization for much of the data collected. All 
array elements operated without failure, and the array 
remained operational throughout its deployment. 

A new mooring and tether design gave the array 
increased Stability. The mean measured array tilt was three 
tenths of a degree. This is meaningful considering that 
previous vertical arrays suffered tilt in excess of ten 
degrees thereby degrading system performance. 

2. Modal Decomposition 

Mode filtering performance is highly dependent on 
array geometry and the normal mode- shapes. Modal 
decomposition requires knowledge of upper and lower boundary 
conditions, sound velocity, and density. The algorithm used 
for modal decomposition is located at the appendix of this 
thesis. (Ref. 13] The decomposition algorithm employs a 
finite differencing routine to approximate the initial eigen 
value problem. Implicit to the routine is the assumptions of 
a pressure release surface and a perfectly rigid lower 
boundary. The bottom is modeled using constant sound speed 
and density characteristics. There is no limitation to the 
quantity or complexity of bottoms except for processing time, 
which increases as the square of the number of points in the 


profile. 


Z.5 


Sound and density profiles used to calculate 
eigenvalues and shapes for the Barents Sea were determined 
from CTD data, and estimations of bottom sound speed and 
density. All processing is based on a single bottom model 
with a constant sound speed of 1667 m/s, and density of 1.83 
g/cm’, taken from Clay and Medwin [Ref. 14]. 

Figure 6 is the Barents Sea sound velocity profile at 
the location of the array. The sound velocity indicates the 
existence of a well defined mixed layer followed by a strong 
negative sound speed gradient extending to the bottom. Due to 
the negative sound speed gradient, the majority of acoustic 
energy is expected to be trapped well below the surface with 
significant bottom interaction. 

3. Normal Modes at the Vertical Array 

Figures 7 to 12 are the mode shapes for the Barents 
Sea array at the 224 Hz source center frequency. Modes one 
through four have turning depths well beneath the surface and 
are expected to contain the majority of acoustic energy based 
on the strong negative sound speed gradient. Modes five 
through ten also turn beneath the surface, but should have 
less beamformed energy than the first four modes because there 
is significant modal energy outside the array’s aperture. 
Modes ten through twenty interact with the surface. These 
modes should have much less energy due to surface scattering 


effects. 
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Figure 7 Sound speed profile at the receiving array. 


Bottom sediment starts at 275 meters modeled by constant op 
1.83 g/cm’, and sound speed 1677 m/s; The bottom is assumed 


rigid below 420 meters. 
4. Modal Beam Pattern 
There are many measurements of performance for 
Classical arrays. Plane wave beamforming is a well understood 


orocess, and performance of a plane wave beamformer can be 
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Figure 9 Normal modes 5 through 8. 
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Figure 10 Normal modes 9 through 12. 
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Figure 11 Normal modes 13 through 16. 
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Figure 12 Normal modes 17 through 20. 
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easily determined through a multitude of available measures 
including: theoretical array gain, estimated array gain, and 
the steering beam pattern. Array gain experienced by modal 
beamforming cannot exceed the theoretical gain for a classical 
beamformer. However, there are additional issues to address 
concerning tomographic resolvability for mode arrivals. 

One such factor, is the cross talk, or the recognition 
differential between modes. To quantify this affect, 
synthetic signals modeled from individual modes were 
beamformed. The modal Signals are created using the frequency 
domain eigen solutions at each hydrophone depth. To duplicate 
the coded signals frequency structure, a Blackman Window is 
applied to the modal signal’s frequency components. The 
square of the beamformed signal is proportional intensity of 
the mode. The peak values of the beamformed modes are squared 
and normalizing by the corresponding modal signal, forming the 
modal recognition kernel. The transpose of the recognition 
kernel gives the modal beam pattern. 

Figures 13 to 18 are the modal beam patterns specific 
to the Barents Sea array geometry and ocean structure. The 
following is a summary of my conclusions and conjectures on 
the Barents Sea array modal beam pattern: 

® Modes one through three are well resolved. 


® Modes four and above are not spatially well resolved due 
to the size of the array’s aperture. 


oS 


Above mode three resolvability must be handled on a case 
to case bases, taking into account the arrival structure 
and relative intensity of each mode. 


Assuming there iS no appreciable modal energy above mode 
twenty, Spatial aliasing effects are non-existent for this 
array. 


Above mode four, there is a recurring beam pattern which 
indicates that the closest modes are the hardest for the 
beamformer to differentiate. Mode four beam pattern shows 
the transition from well resolved to less resolvable 
arrivals. 


Increasing array length to 50 meters below the surface 
would increase the coverage of the array’s aperture and 
guarantee that at least the first six modes are spatially 
resolved. This can be accomplished through increasing the 
number of hydrophones, or by expanding element spacing, 
subject to spatial aliasing. 
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VI. EVALUATION OF BARENTS SEA DATA 


A. DECODING OF RAW ACOUSTIC DATA 





The participants of the Barents Sea Polar Front Experiment 
performed some near-real-time signal processing during the 
experiment. Some of the signal processing included decoding 
of the 224 Hz source. Figure 18 is the decoded signal taken 


from the bottom hydrophone of the array. 


Acoustic time series after decoding 


magnitude squared 


80 


tume (sec) 





Figure 18 Acoustic time series after BM decoding for the 
bottom hydrophone in the array [Ref. 7]. 


The BM sequence removal algorithm performed very well. 


Theoretical decoding gain is based on the number of digits in 
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the coded signal. The 224 Hz source signal has a theoretical 
decoding compression gain of 18 dB which compares favorably to 
initial estimates. Signal gain estimates where determined by, 


SNR 
GAIN = 10109 (Sean , (6.1) 
MINE sy cyge 


Meet. 15]. 

Coherent averaging is a useful method for determining the 
acoustic arrival structure across the array. Figure 19 is the 
coherently averaged signal across the array, post m-sequence 
removal . The coherently averaged signal indicates the 
presence of at least two separate arrivals. This effect is 
most likely do to the time spreading of the signal between the 
propagation paths. The first arrival has the bulk of acoustic 
energy and a complex structure. Clearly this energy is formed 
by the bottom bounce propagation path. The second arrival 
lacks the energy and fine definition of the first caused by 


Fncreased bottom and surface interaction. 


B. RESULTS FROM THE BROADBAND MODAL BEAMFORMER 

The results from broadband modal beamforming of the 
Barents Sea data are quite impressive. The mode arrival 
shapes are very distinctive, which simplifies arrival time 
estimation. All times for Figures 21-27 are relative to the 
record tape, and do not represent the travel from source 


emission. Figure 20 shows the mode one arrival structure for 
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Figure 19 Coherent averaging after BM decoding for the 
WHOI-NPS array [Ref. 7]. 


the first four sequences of decoded pulses. Figures 21 
through 23 are the first five modes pertaining to the first 
four sequence arrivals of Figure 18. The first five modes are 
the strongest arrivals. This is expected because the majority 
of modal energy is within the array aperture. Mode three is 
the strongest arrival indicating the majority of acoustic 
energy is in the lower 125 meters of the water column. Mode 
two is particularly interesting due to its multiple arrivals. 
This pattern might be attributed to coupling caused by the 
sharp change in sound speed gradient at the front. The mode 
two arrival is both the smallest in total energy and the 


earliest arrival of the first five beamformed modes. 
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Figure 20 Mode one, beamformed and decoded arrival 
SeEructure. 


Modes six through ten, Figure 25, appear after the initial 
arrival structure. These modes contain less energy and have 
later arrival times. A large amount of modal energy exists 
outside the array aperture and is at least a partial factor to 
the decrease in beamformed energy. 

Modes eleven through twenty, Figures 26 and 27, show a 
sharp drop off in energy levels which cannot be completely 
contributed to the array aperture. This sudden decrease of 
beamformed energy 1S most likely caused by high transmission 
Hess Caused by surface scattering. It is interesting to note 


that arrivals for modes sixteen and seventeen indicate 
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increased energy levels. This is a curious result that 
presently eludes explanation. 

Broadband modal beamforming of the BSPFEX data provided 
extremely clear arrival structure that could easily be 
mistaken for modeled data. Part of the reason that the 
results are this good is the extremely high signal to noise 
ratios and the low propagation loss. Clearly, the first ten 
modes are the most Significant to acoustic tomography. Mode 
resolvability is dependent upon the mode’s arrival time, the 
modal beam pattern and spatial aliasing. Resolvability issues 
have not yet been completely sorted out but may indicate that 
only the first few modes are reliable enough to use in the 


tomographic inverse. 
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Figure 21 First arrival sequence for modes one through five. 
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Figure 22 Second arrival sequence for modes one through five. 
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Figure 24 Fourth arrival sequence for modes one through five. 
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Arrival structure for modes six through ten. 
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Figure 26 Arrival structure for modes eleven through fifteen. 
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VII. SUMMARY AND FUTURE WORK 


A. CONCLUSIONS AND SUMMARY 

Preliminary results indicate that broadband modal 
beamforming is a useful tool in shallow water acoustic 
tomography. Ray acoustic tomography in shallow water is 
complicated because the earliest ray arrivals, containing the 
bulk of the acoustic energy, are unresolved. Unresolved ray 
paths are well described by the lower mode arrivals of the 
beamformed array. Broadband modal beamforming produces 
distinct mode energy arrivals simplifying estimation of 
arrival times. 

Improvements made by WHOI and NPS to the vertical array 
design proved very successful. The Dual Benthos interrogators 
provided a more reliable source for measuring array geometry 
than previous methods. Array stability was vastly improved by 
its new moored design and should be used as a benchmark for 
the design of future arrays. The measured array tilt during 
the experiment was a nominal three tenths of a degree and 
minimized the correction required for array tilt. 

Based on recommendations from Crocker [Ref. 4], a fast 
third order interpolation routine was implemented to reduce 
processing time. The new interpolator was derived using the 


Taylor Approximation for a third order polynomial [Ref. 16]. 


a2 


This algorithm is approximately five times faster than 
previous interpolator based on Neville’s algorithm 


[Ref. 17]. 


B. RECOMMENDATIONS 

The broadband modal beamformer requires each mode to be 
processed individually. Most of the time incurred from 
processing is due to the Fast Fourier Transform algorithm. 
Processing all modes simultaneously would reduce the number of 
Fast Fourier Transforms by the number of desired modes. 
Significant Savings in processing can be realized by adopting 
this technique. Calculating the modal beam pattern was an 
additional time sink. Since each mode pattern must be 
processed for every mode, processing time goes as the square 
of the number of modes desired in the beam pattern. IMeng 
twenty modes, the beamformer must be run four hundred times to 
get the beam pattern for just one sound velocity profile. If 
all the modes where processed simultaneously, the modal beam 
pattern could be an added option to the program reducing the 
processing required by the user. 

Decoding of the m-sequence coded signal occurs after 
beamforming. The decoding process could have been 
accomplished prior to beamforming but would have increased 
processing time by the number of elements in the array. The 
output from the m-sequence removal algorithm is complex 


numbers thereby doubling the amount of memory required to 


oo 


Store data. Despite this, it may still be desirable to decode 
the data first so that plane wave and modal beamforming can be 
accomplished using the same data base. To perform sequence 
removal prior to beamforming, the decoder algorithm must be 
adapted to handle multi-channel data, and an interpolation 
routine should be added to resample data prior m-sequence 
removal. The modal beamformer will also need to be updated to 
handle complex input and to perform mode filtering on 
demodulated signals. 

Design of the next vertical array should incorporate 
estimations of the modal beam pattern. Modeling of the 
acoustic environment will give an estimate of the modal beam 
pattern. Using this information, the number of array elements 
and geometry can be optimized to achieve the desired 


performance level. 
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APPENDIX 


A. BROADBAND MODAL BEAMFORMER ALGORITHM 
The broadband modal beamformer was designed for a robust 
environment. Program portability allows it to be used on any 
computer system supporting C or quick C software. The program 
is flexible enough to cope with most vertical line array 
designs. A large number of user parameters are incorporated 
in the computer architecture to enhance portability. The 
following is a summery of user parameters employed. 
1. Program Parameters 
@® UNIX VERSION: a compiler flag used to determine the 
operating system, either UNIX or ANSI. 
® ASCII: ASCII and BINARY options determine the beamformed 
output type. These options are inter-linked, and never 


can be activated at the same time. 


® SIGNAL: enables time series beamformed output. This 
parameter should be left "ON". 


® LOWER _ SENSOR: enables lower sensor tilt data if available. 


® VALIDATE: this option flags the program to output several 
validation files. If selected, the program will output 
files containing array geometry, steering delays, 
hydrophone weights at the lowest frequency in the band, 
and estimated mode speed. 


@® ERROR_ESTIMATE: flags to produce the estimated error 
caused by interpolation. The error estimate is stored in 
a file specified by the user. 


® VERBOSE: if selected, this option will cause the program 
to provide a status report during operation. Information 
is sent to standard output. This option is useful for 
debugging and for short data sets. 


@® ON/OFF: logical switches for program control. 
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INTERPOLATE: specifies which interpolation routine to use. 
The available interpolators are fastpoly3, polint and 
Bat dite Fastpoly3 is the fastest of the three 
interpolators. Other interpolators can be added by the 
user if desired. 


ORDER: determines the order for interpolation. Polint and 
ratint are multi order algorithms, while fastpoly3 can 
only perform third order interpolation. If the order of 
interpolation is not three and the fastpoly3 interpolator 
is selected, the program will terminate on a standard 
error. 


STEP: indicates the number of points on either side of a 
sequence for derivative estimates. 


TINY: smallest number allowed for floating point 
operations. 


PI: trigonometric constant. 
RADIAN: degree to radian conversion constant. 


OFFSET: correction term for the difference between tilt 
sensor depth and the depth of the first hydrophone. 


DELTA R: array element spacang an meters. 
CTD OFFSET: used to correct for offset Caused by she ems 
measurement. This offset is reflected in the eigen 


solutions. 


SAMPLE DELAY: corrects for sampling skew caused by 
sampling data from a single data bus (ms). 


SSP LENGTH: dimensioning term based on the maximum number 
of points in a single eigen function. 


EIGVAL LENGTH: dimensioning term based on the maximum 
number of frequencies contain in the eigen value file. 


LOOK DIRECTION: desired direction for beamforming in 
degrees true. 


TILT_BUFFER: dimensioning factor based on the maximum 
number of tilt observations in the tilt input file. 


BUFFER_TIME: dimensioning term. Must be an integer and 
greater then the FFT TIME + 2. 


F SAMPLE: the integer sampling frequency. 


ae: 


———————————<_—_——— Ss 


@® CHANNELS: total number of array elements. 


@ F CARRIER: 


the desired carrier frequency to beamform. 


Note that this should match the center frequency of the 
eigen function and value input files. 


Seer. TIME: .amount of data, 


FFT. This should be the sequence length of the desired 


beamformed signal. 


@® FFT LENGTH: the radix two FFT size. 


in seconds, 


be larger or equal to FFT _TIME*F SAMPLE. 


® SWAP: macro used by the FFT algorithm. 


Zs 


[xxx 


+ + + + € FF F HF F HF HH F KF KH OF 


* 
* ke / 


#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 


Outputs a data file. 


Beamformer Source Code 


PROGRAM: BEAMFORMER vsn 5.0 
NAME : 
WRITTEN BY: Glenn A. Omans II / Steven Crocker 


BBbeam.c 


LAST UPDATE: September 28, 


This program takes input from various data files and the user. 
The inputs are a number of channels of digital 
acoustic data, and information regarding the physical characteristics 
and geometry of the receiving array. 
in the form of normal mode eigenfunctions and eigenvalues at the 
receiving array are required to operate this beamformer. 
is a single channel of acoustic data. 


UNIX VERSION 

ASCII ON 

BINARY OFF 

SIGNAL ON 
LOWER_SENSOR OFF 
VALIDATE OFF 
ERROR_ESTIMATE OFF 
VERBOSE OFF 

ON 1 

OFF 0 

INTERPOLATE fastpoly3 
ORDER 3 

STEP 1 

ten, 2..0e-25 

Pie 3, 14159265359 
RADIAN 57.2957795131 
OFFSET 0.0 

DELTA _R 10.0 


PURPOSE: Broadband Modal Beamforming of a Vertical Array. 


Vo 22 


jie 


allotted to each 


The FFT LENGTH must 


Additionally, environmental data 


The output 


/* either ANSI or UNIX * / 
/* select output mode ON or OFF * / 
/* select output mode ON or OFF * f 
/* either ON or OFF * / 
/* either ON or OFF i 
/* either ON or OFF aye 
/* either ON or OFF =) 
/* either ON or OFF wif 
/* logical "switch" z/ 
/* logical "switch" = 
/* either polint, fastpoly3 or ratint */ 
/* order of interpolator (odd) */ 
/* number of steps for derivatives ay 
/* prevents division by zero ig 
/* for freq to omega conversions as 
/* for degree to radian conversions * / 
/* dist btwn upper sensor and phone #1 */ 
/* array element spacing = / 


=f 


#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 


#include<stdio.h> 
#include<malloc.h> 
#include<math.h> 


CTD_OFFSET 0.0 /* aiff btwn ctd depth inc & 1st depth © 





SAMPLE DELAY 0.003 /* array sampling delay (ms) 
SSP_LENGTH 2500 /* max number of pts in eigunfunction 7 
EIGVAL_ LENGTH 230 /* max number of eigenvalues 1 
LOOK_DIRECTION 270.0 /* direction from which signal arrives + 
TILT BUFFER 120 /* max length of tilt data vectors 7 
BUFFER_TIME 7 /* input buffer length in seconds (int)? 
F SAMPLE 1600 /* sampling frequency (int) ? 
CHANNELS 16 /* number of channels processed * 
F CARRIER 224.0 /* carrier frequency 
FFT_TIME 3.9375 /* time alloted to each fft * 
FFT LENGTH 8192 /* radix 2 <= (FFT_TIME-2)*F SAMPLE * 


SWAP (a,b) tempr=(a) ; (a) =(b) ; (b) =tempr 


#if defined ( ANSI ) 
#include<float.h> 
#include<stdlib.h> 
int getInput (void) ; 
int putOutput (void) ; 
int processTilt (float **x, float **y, float **z, float *delayClock) ; 
int processModes (float **z, float **weight, float *ptrcC) ; 
int dydx(float *x, float *y, float *ddx, int points) ; 
int fastpoly2 (float *xa, float *ya, int n, float x, float *y, float *dy) . 
int fastpoly3 (float *xa, float *ya, int n, float x, float *y, float *dy) , 
int polint (float *xa, float *ya, int n, float x, float *y, float *dy); 
int ratint (float *xa, float *ya, int n, float x, float *y, float *dy) ; 
int realft (float *data, int n, int isign); 
int fourl(float *data, int nn, int isign); 
int window(float *data, int N); 
float *vector(int length) ; 
float **matrix(int row, int col); 
int free_matrix(float **m, int row) ; 
void ExitOnError (char error txt[]) ; 
#Helif defined ( UNIX ) 
int vector (), 
matrix(), 
getInput (), 
putOutput (), 
processTilt(), 
processModes (), 
fastpoly2(), 
fastpoly3(), 
poling), 
ratint(), 
dumpSpectrum(), 
realft(), 
fourl () t 
free _matrix(), 
ExitOnError () ; 


#tendif 


/* Global Variable Declarations */ 
float **inSOUND, *outSOUND, *MeanSqError, Max=0.0, Min=1.0e25; 
float fmax, fmin; 
int LastDelay, Minute=1, firstBuffer=1; 
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FILE *fpinSound, *fpOutSound, *fpOutSpectrum, *fpMSE; 


main () 


/* 


char fileDelay[80], fileModes[80], fileArray [80] ; 
FILE *fpV1; 
mitt, J, Kk, @, Ltime, touff£, kount, *shift; 
float **x, **y, **z, *indx, *samples, arg, ans, err, *delay, 
*delta, **weight, *pwrSpectrum, Cgroup, *Xf, 
*sumXf, **buffSOUND, df, freq, *delayClock; 
double Time; 
/* Memory Allocation and Tilt Data Processing */ 
x= (float**)matrix (CHANNELS, TILT BUFFER) ; 
y= (float**) matrix (CHANNELS, TILT_BUFFER) ; 
z=(float**)matrix (CHANNELS, TILT_BUFFER) ; 
delayClock=s (float*) vector (TILT_BUFFER) ; 
processTilt(x, y, z, delayClock) ; 


printf ("\nDelay Clock Times\n\n") ; 
for (i=0; i<=LastDelay; i++) printf ("%i\t%g\n",i,delayClock[i]); */ 


if (VALIDATE==ON) 
printf ("\nEnter file name for array validation data: "); 
if ((scanf ("%s",fileArray) ) ==EOF) 
ExitOnError ("Fatal error in scanf()"); 


fpVl=fopen(fileArray,"wt") ; 
if (fpVl==NULL) ExitOnError ("Error opening validation file") ; 


ponte rt (epvl ae ehannel \t yt X\t\toy \t\t 2\n")-; 
for (i=1;i<=LastDelay;i++) 


gexrintt (fpVil;, MMINUTE: Si\n", 1); 
for (j=1;3<=CHANNELS ; j++) 


fp eens £ {fp V1; Csr ee Sf \ € S\N e 
eae x3) £1) ,y 5) lil,z zi [1]); 
[= EOr = / 
fclose (fpV1) ; 
} /* if */ 


/* Memory Allocation*/ 

weight=(float**) matrix (CHANNELS, FFT_LENGTH) ; 

indx= (float*) vector (ORDER+1) ; 

delay= (£loat*) vector (CHANNELS) ; 

delta=(float*) vector (CHANNELS) ; 

outSOUND= (float*) vector ( (BUFFER | TIME-1)*F SAMPLE) ; 

Xf = (float*) vector ( (FFT_LENGTH+1) *2) ; 

sumxf= (float*) vector ( (FFT_LENGTH+1) *2) ; 

buf f SOUND= (float**) matrix (CHANNELS, (BUFFER_ TIME -1) *F_SAMPLE+1) ; 

inSOUND= (float**)matrix (CHANNELS, BUFFER TIME*F SAMPLE) ; 

ae ( (shaiete= (int*) malloc ( (CHANNELS+1) *sizeof (int)) )= =NULL) 
ExitOnError ("Memory allocation failure for shift[]. NN ag 


Se) 


for (i=1; i<=ORDER+1; i++) indx[i] =(float)i; 
cs 2 ( ER R O R._ ES T IOMOA TT 2 =e] 70 eae 
MeanSqError= (float*) vector ( (BUFFER_TIME-1) *F_SAMPLE) ; 


/* Mode Data Processing */ 
processModes(z, weight, &Cgroup) ; 


/* Calculate delays, shifts, ect */ 
for (1=1; 1<=CHANNELS; i++) 


delay [i] =x[i] [Minute] /Cgroup + /* Time delay * / 
(float) (i-1)*SAMPLE DELAY/1000.;  /* Sampling delay */ 
shift [i] =(int) (delay [i] * (float) F_ SAMPLE) ; /* # of samples */ 


/* fraction of 1 sample 
oll 
delta [1] =delay [i] * (float) F_SAMPLE- (float) shift [1] ; 
\ /* fora 


while (Minute<=LastDelay) 


if (VALIDATE==ON) 
if (firstBuffer) 
printf ("\nEnter file name for delay validation: ")j; 


if ((scanf ("t$s",fileDelay) ) ==EOF) 
ExitOnError ("Fatal error in scanf()"); 


if ((fpVl=fopen(fileDelay, "wt")) == NULL) 
ExitOnError("Error opening validation file."); 
ey ol alge 7) 
else 
if ((fpVl=fopen(fileDelay, "at")) == NULL) 


ExitOnError ("Error opening validation file."); 
} /* else */ 


Eprintt (fpV1, "MINUTE: 62 ne Maines 


fprinee (epvi, 
"Channel\t delay\t\t int ehift\t traction cl I shitt ne 


for (i=1;i<=CHANNELS; i++) 


fprintf(fpV1, "%i\t te\t %i\té e\n", 
1,delay [1] ,shift[i] ,delta[i]) ; 
/* for */ 
fclose (fpvV1) ; 


if (firstBuffer) 


printf ("\nEnter file name for phone weights and group speed: 
i) 
if ((scanf ("ts", £ileModes) ) ==EOF) 
ExitOnError ("Fatal error in scanf()"); 


cafe ( (£pV1=fopen (fileModes, "wt")) == NULL) 
ExitOnError ("Error opening validation file.\n") ; 
fprintf (fpvl, "Group speed £On F_CARRIER PS = 


tq\n\n\nh"; Cqroup)y, 
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Pe tt Lixreteurter */ 
else 


if ((fpVl=fopen(fileModes,"at")) == NULL) 
ExitOnError ("Error opening validation file.\n") ; 
fprintf (fpVv1,"\nHydrophone weights for minute %i\n",Minute) ; 
} /* else if not firstBuffer */ 
fprintf(fpV1, "Channel \tWeight\n") ; 


for (i=1;1<=CHANNELS;i++) 
Eprintce (fpVi, %¥i\tse\n",i,weightiy) [1 lee 


fclose (fpV1) ; 


} /* iff VALIDATE */ 


Time = 1; /* First second of time is 


chopped */ 


delay 


for (n=1; n<=60/(BUFFER_TIME-2); n++) /* for one Minute */ 


getInput () ; 


if (firstBuffer) /* Produce first output buffer */ 


{ 
if (VERBOSE) printf("Interpolation first Buffer\n") ; 
for (i=1; i<=F_SAMPLE; i++) 
OUutSOUND [1] =0.0; 
if (ERROR_ESTIMATE==ON) MeanSqgError [i] =0.0; 
lie *ebor a) 
for (1=F SAMPLE+1; 1<=(FFT_TIME+1) *F_SAMPLE; 1++) 


outSOUND [1] =0.0; 
if (ERROR_ESTIMATE==ON) MeanSqError[i]=0.0; 


Time+=(1/ (float) F_SAMPLE) ; 
if (Time >= delayClock [Minute] ) 
Minute++; 


/* Mode Data Processing */ 
processModes(z, weight, &Cgroup) ; 


/* Calculate delays, shifts, ect */ 
for (k=1; k<=CHANNELS; k++) 


ole 


Sampling delay */ 


samples */ 


delay [k] =x[k] [Minute] /Cgroup + ) tame 
(float) (k-1)*SAMPLE DELAY/1000. ; /% 
shift [k] =(int) (delay [k] * (float) F_SAMPLE) ; /* # of 
/* 
fraction of 1 sample */ 
delta [k] =delay [k] * (float) F_SAMPLE- (float) shift [k] ; 
} /* for */ 


} /* if (Need new shift delays) */ 


ae (j=1; j<=CHANNELS; j++) 


Syl 


if (delay [j] >=0.0) 


arg=(float) (ORDER+1) /2.0+(1-delta[j]); 

samples = &inSOUND [Jj] [i-shift [j] - (ORDER+1) /2] ; 
} /* if */ 
else if (delay[j] <0.0) 


arg=(float) (ORDER+1) /2.0+deltal[j] ; 
samples = &inSOUND[j] [i-shift [j] - (ORDER+1) /2+1] ; 
} /* else if */ 


INTERPOLATE (indx, samples, ORDER+1, arg, &ans, &err) ; 
buffSOUND [j] [i-F_SAMPLE] =ans; 


if (ERROR_ESTIMATE==ON) 
MeanSqgError [i] =MeanSqError [i] +err*err; 
ay Oi 
if (ERROR_ESTIMATE==ON) 
MeanSqError [i] =MeanSqError [i] / (float) CHANNELS ; 


yey Ora 
ES PGT erty 
else /* Produce subsequent output buffers */ 


ae (VERBOSE) printf ("Interpolating subsequent output 
buffers \n' 


for (i1=F_SAMPLE+1; i<=(FFT_TIME+1) *F_ SAMPLE; i++) 


outSOUND [i-F SAMPLE] =0.0; 
if (ERROR_ESTIMATE==ON) MeanSqError[i-F_ SAMPLE] =0.0; 


Time+=(1/(float)F SAMPLE); /* time processed with first 
second chopped */ 


if (Time >= delayClock [Minute] ) 
Minute++; 


/* Mode Data Processing */ 
processModes(z, weight, &Cgroup) ; 


/* Calculate delays, shifts, ect */ 
for (k=1; kK<=CHANNELS; k++) 


delay [k] =x[k] [Minute] /Cgroup + /* Time 
delay */ 
(float) (k-1)*SAMPLE DELAY/1000. ; ihe 
Sampling delay */ 
shift [k] =(int) (delay [k] * (float) F_SAMPLE) ; /* #¥og 
samples */ 
/* 


fraction of 1 sample */ 
delta [k] =delay [k] * (float) F_SAMPLE- (float) shift [k] ; 


} /* for */ 
} /* if (Need new shift delays) +*/ 
for (j=1; j<=CHANNELS; j++) 

if (delay [j] >=0.0) 


arg=(float) (ORDER+1) /2.0+(1-delta[j]); 
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samples = &inSOUND[j] [i-shift[j] - (ORDER+1) /2] ; 
ax iE */ 
else if (delay[j] <0.0) 


arg=(float) (ORDER+1) /2.0+delta[j] ; 
samples = &inSOUND [Jj] [i-shift [j] - (ORDER+1) /2+1] ; 
/* else if */ 
INTERPOLATE (indx, samples, ORDER+1, arg, &ans, &err); 
buffSOUND [j] [1-F_SAMPLE] =ans; 
if (ERROR_ESTIMATE==ON) 


MeanSqError [i1-F SAMPLE] =MeanSqError [1-F_SAMPLE] +err*err; 
J J* for =) 
if (ERROR_ESTIMATE==ON) 





MeanSqError [1-F_SAMPLE] =MeanSqError [1-F_SAMPLE] / (float) CHANNELS ; 
| yf for */ 
: } /* else */ 


/* take fft multiply by weights & sum in frequency domain */ 
df= (float) F_SAMPLE/ (float) FFT LENGTH; 
if (VERBOSE) printf("\nTaking FFT\n") ; 
Poni i=t,  tec=2*FPT LENGTH; i++) sumXf[i]=0; 
for (j=1;j3<=CHANNELS ; j++) 
for (1=1;1<=FFT LENGTH; i++) 
if(i <= FFT_TIME*F_ SAMPLE) 


Xf (2*i-1] =buffSOUND [j] [1] ; 
Xf [2*i] =0; 


else 


Xf [2*1i-1] =0; 
Xf [2*1i] =0; 


Poh Silopa 3 a7) 
fourl (Xf, FFT_LENGTH, 1) ; 
k=1; 


for(i=l1; i <=FFT_LENGTH-1; 1i+=2) 


freq=df* (float) (i-1)/2; 
if (freq>=fmin && freq<=fmax) 


/* positive frequencies */ 

sumXf [1] =sumXf [i] +Xf [i] *weight [Jj] [k] ; /* real part */ 
sumXf [1+1] =sumXf [i+1] +Xf [i+1] *weight [5] [k]; /* imag part */ 
/* negative frequencies */ 


sumXf [2*FFT_LENGTH-i] =sumXf [2*FFT_LENGTH-i] +Xf (2*FFT LENGTH-i] *weight [j] 
[k] ; 7 


sumXf [2*FFT_LENGTH-1+1] =sumXf [2*FFT_LENGTH-i+1]+Xf[2*FFT LENGTH-i+1] *wei 
oie) ] (kl; = 
k++; 
/* if (frequency within Bandwidth) */ 
he = SE Or ee 
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} o/* feria, 
if (VERBOSE) printf("\nTaking invFFT\n") ; 
fourl (sumxf , FFT_LENGTH, -1) ; 


if (firstBuffer) j=F_SAMPLE+1; 
else j = 1; 
for(i=1; 1 <=F_SAMPLE*FFT_TIME;1i++) 


outSOUND [53] =(1/ (float) FFT_LENGTH) *sumXf [2*i-1] ; 





if (fabs (outSOUND [j])<Min) Min=fabs (outSOUND [j]) ; 
if (fabs (outSOUND [j]) >Max) Max=fabs (outSOUND [3]) ; 
J++; 


) o/* forms, 


/* clear summer for next set of data */ 
for(i=1; 1 <=2* FFT LENGTH; 1+) se sumxt (2) =C; 


if (VERBOSE) printf("\nWriting Data to file\n") ; 
if (SIGNAL==ON) putOutput () ; 
firstBuffer=0; /* set firstBuffer to false */ 

| eons), 


if (VERBOSE==ON) 
printf ("\t%i minutes of input data processed.\n", Minute) ; 


if (VERBOSE==OFF && Minute==1) 


printf ("\n\nAll user input has been accepted.\n") ; 

printf ("Program is in SILENT mode.\n") ; 

printf("If desired, put this job in background now.\n") ; 
Weir + 


} /* while not next minute*/ 


if (VERBOSE==ON) 


printf ("EXECUTION COMPLETE: End of tilt data encountered\n") ; 
printf ("Maximum magnitude encountered was: %te\n",Max) ; 
printf ("Minumum magnitude encountered was: %e\n",Min) ; 


fclose (fpInSound) ; 

if (SIGNAL==ON) fclose (fpOutSound) ; 

if (ERROR_ESTIMATE==ON) fclose(fpMSE) ; 

exit (0) ; 
[ee HEKKRER KEKE IKKE KREERKEEKEEKKERERKKKKEKKKKKKKKKKKK KKK K KK KKK KKK / 
[RE RERREAERER EEK EA RES © = Se NE) main tok kok kk kk kkk kkk kk kk ek  / 
[BREE HKK KEENER EREEEREEEEREREREREKRE RE EK RE RIK Re eee ee 

[xxx 


* FUNCTION: getInput () 

* 

* This function handles all acoustic input. It also provides one of 

* two normal process terminations available in the program. (The other 
* is located in main().) 

* 

* Arguments: none 

* 

* Return value: 0 

* 
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Be Funceions called: vector () Ex1itOnError () 
* 

e Definitions called: ANSI UNIX 

‘g F SAMPLE CHANNELS 

‘a BUFFER_TIME VERBOSE 

a ERROR_ESTIMATE SIGNAL 

* 

* Global variables called: 1nSOUND [] [] Min 

- firstBuffer Max 

a fpInSound fpOutSound 
a fpOutSpectrum fpMSE 

* 

* Significant memory allocation: diskBuffer [] 

* 

kkkkk 


#if defined ( ANSI ) 
int getInput (void) 
#Helif defined ( UNIX ) 
getInput () 

#Hendif 


int i, j, buffer, items; 
float *diskBuffer; 
char fileName [80]; 


if (firstBuffer) 
printf ("Enter file name for input acoustic data: "); 


if ((scanf("%s", fileName) ) ==EOF) 
ExitOnError("Fatal error in scanf()"); 


preenet«'' \a\n") ; 


if ((fpInSound=fopen(fileName, "rb")) == NULL) 
ExitOnError ("Error opening INPUT ACOUSTIC data file."); 


eet */ 


/* Memory Allocation */ 

if (firstBuffer) buffer=(2+FFT_TIME) *F SAMPLE*CHANNELS ; 
else buffer=(FFT_TIME) *F_SAMPLE*CHANNELS; 
diskBuffer= (float*) vector (buffer) ; 


items=fread((char*) (diskBuffer+1), sizeof (float), buffer, fpInSound) ; 
if (VERBOSE) printf ("items = %i\n",items) ; 
if (items==buffer) /*continue*/; 


else if (ferror(fpInSound) != 0) 
ExitOnError ("Error encountered while reading input acoustic data") ; 
else if (feof (fpInSound) != 0) 


if (VERBOSE==ON) 


Brimet (Vac. ee ek Oe AREER RARER EKER ERK KKKKKKKKKKKAKEKKEEK\ YN) ; 


jongueenene | aAN End of File reached: EXECUTION COMPLETE\n") ; 
printf ("\t\t%i minutes of data processed\n",Minute-1) ; 

printf ("\t\t%i bytes of data discarded\n", items*sizeof (float) ) ; 
printf ("\t\tMaximum magnitude encountered was: %e\n",Max) ; 
printf ("\t\tMinimum magnitude encountered was: ¢¥e\n",Min) ; 
tet Ge End of File reached: EXECUTION COMPLETE\n") ; 


o 
BIN CE (1107 I III III II IO I tok tata INDI); 
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} 


fclose (fpInSound) ; 
if (SIGNAL==ON) fclose(fpOutSound) ; 
if (ERROR_ESTIMATE==ON) fclose(fpMSE) ; 
exit (0); 
} /* else if */ 
else ExitOnError ("Unknown error handling acoustic input file."); 


for (i=1; i<=(2+FFT_TIME)*F SAMPLE; i++) 
for (j=1; j<=CHANNELS; j++) 
if (firstBuffer) 
inSOUND [j] [i] = diskBuffer [CHANNELS* (i-1)+j] ; 
es se 
else 


if (i<=2*F_ SAMPLE) 


inSOUND [j] [1] =inSOUND [j] [1+ (int) ( (float) FFT_TIME* (float) F_SAMPLE) ] ; 


else 
inSOUND [j] [1] =diskBuffer [CHANNELS* (i-2*F_SAMPLE-1) +]j] ; 
} /* else */ 
/* for */ 
} /* for */ 


/* Deallocate Memory */ 
free ((char*) diskBuffer) ; 
return( 0 ); 


[RRR KKK RRR KR KR RK RK / 


[/***** END getInput wk ke / 


[RRR KK RK KR RK KK RK CK / 


[eee 

* FUNCTION: putOutput () 

* 

* This function handles all acoustic output. Additionally, it outputs 
* the estimated mean squared error from the interpolators (if enabled). 
* 

* Arguments: none 

* 

* Return value: 0 

* 

* Functions called: ExitOnError () 

te 

* Definitions called: ANSI UNIX 

* F_SAMPLE BUFFER TIME 
* ERROR_ESTIMATE ASCII 

* BINARY 

* 

* Global variables called: outSOUND [] MeanSqgError [] 
* firstBuffer fpOutSound 

* fpMSE 

* 

* Significant memory allocation: none 

* 

tk ete / 


#if defined ( ANSI ) 
int putOutput (void) 
Helif defined ( UNIX ) 
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putOutput () 
_ 


fete 1, cut 
char fileName[80], *mode; 


me (fsrstBuftfer) 


{ 
if (ASCII==ON) mode="wt"; 
else if (BINARY==ON) mode="wb"; 
Suc=1 ; 
printf ("Enter file name for output data: ") ; 
if ((scanf("%s", fileName) ) ==EOF) 
ExitOnError ("Fatal error in scanf()"); 
printf ("\n\n"); — 
if ((fpOutSound=fopen (fileName, mode)) == NULL) 
ExitOnError ("Error opening OUTPUT data file."); 
if (ERROR_ESTIMATE==ON) 
printf ("\nEnter file name for interpolator error estimate: "); 
if ((scanf ("%s", fileName) ) ==EOF) 
ExitOnError ("Fatal error in scanf") ; 
if ((fpMSE=fopen (fileName, mode)) == NULL) 
ExitOnError ("Error opening error.dat") ; 
Pe ie / 
aac */ 


else cut=2; 
if (ERROR_ESTIMATE==ON) 
if (ASCII==ON) 


for (i=1; i<=(2+FFT_TIME-cut) *F_SAMPLE; i++) 
fprintf (fpMSE, "sf\n", MeanSqError [i]) ; 
} /* af */ 
else if (BINARY==ON) 


if (fwrite ((char*) (MeanSgError+1) ,sizeof (float) , (2+FFT_TIME-cut) *F_SAMPLE, 
fpMSE) == (unsigned) (2+FFT_TIME-cut)*F_ SAMPLE) ; 
else if (ferror(fpMSE) != 0) 
ExitOnError ("Error encountered writing error data") ; 
else ExitOnError ("Unknown error handling error file."); 
} /* else if */ 
} /* if * / 
if (ASCII==ON) 
for (1=1;1<=(2+FFT_TIME-cut) *F_ SAMPLE; i++) 
fprintf (fpOutSound, "%sg\n", outSOUND[i}) ; 


VE Aas aly 
else if (BINARY==ON) 


if (fwrite ((char*) (outSOUND+1) , sizeof (float) , (2+FFT_TIME-cut) *F_ SAMPLE, 


oy 


+ £¢ + ¢ + + £ + € FE FF FF FF OF FF FF FF HF OH FE HF EN 


fpOutSound) == (unsigned) (2+FFT_TIME-cut)*F SAMPLE) ; 
else if (ferror (fpOutSound) != 0) 
ExitOnError ("Error encountered writing output acoustic data") ; 
else ExitOnError ("Unknown error handling acoustic output file.") ; 
} /* else if */ 
return( 0 ); 


[wR a Kk a ke kk tk tek kok tok ee / 

/***** END putOutput *****/ 

[#8 ee ik te toe doe koe te kde ke koe tok ee / 
kkkkk 


FUNCTION: processTilt () 


This function handles all array tilt data. It calculates the X, Y, Z 
coordinates of each hydrophone as a function of time. The coordinate 
system is oriented such that X points toward the signal "origin" and 
Z points down. 


Arguments: ame y(J (0) 
z() () 
Return value: 0 
Functions called: ExitOnError() matrix() 
vector () free _matrix() 
Definitions called: ANSI UNIX 
DELTA _R LOWER_SENSOR 
RADIAN CHANNELS 


LOOK_DIRECTION OFFSET 
TILT BUFFER 


Global variables called: LastDelay 
Significant memory allocation: xx [J] [] yy [] [] 
zz[) [J eek sie |) 
angle [] udepth [] 
ldepth [] delayClock [] 
* 
we dk / 


#if defined ( ANSI ) 

int processTilt (float **x, float **y, float **z, float *delayClock) 
#telif defined ( UNIX ) 

processTilt (x,y,z, delayClock) 

float **x, **y, **z, *delayClock; 

#tendif 


int sis ai) NOCEOrs 

float *tilt, *angle, *udepth, *ldepth, **xx, **yy, **zz, theta; 
char fileName [80] ; 

FILE *ipl, =f£o2, 


/* Open Data Files */ 
printf ("Enter file name for upper tilt sensor data: "); 


if ((scanf("%s", fileName) ) ==EOF) 
ExitOnError ("Fatal error in scanf()") ; 


eparberera i ere) 


if ((fpl=fopen(fileName, "rt")) == NULL) 
ExitOnError ("Error opening UPPER TILT data file."); 
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if (LOWER_SENSOR==ON) 
printf ("Enter file name for lower tilt sensor data: "); 


if ((scanf ("%s", fileName) ) ==EOF) 
ExitOnError ("Fatal error in scanf()") ; 


PrintLy can \n")., 


if ((fp2=fopen(fileName, "rt")) == NULL) 
ExitOnError ("Error opening LOWER TILT data file."); 
ldepth= (float*) vector (TILT_BUFFER) ; 
} /* if */ 


/* Memory Allocation */ 
Gult—(ricat*) vector (TILT_BUFFER) ; 
angle=(float*) vector (TILT_BUFFER) ; 

udepth= (float*) vector (TILT_BUFFER) ; 

xx= (float**)matrix (CHANNELS, TILT_BUFFER) ; 
yy= (float**)matrix(CHANNELS, TILT_BUFFER) ; 
zz=(float**)matrix(CHANNELS, TILT _BUFFER) ; 


1: /* read upper tilt data */ 
MOEROF=1;> 
while (notEOPF) 
if (fscanf (fp1,"%g tg tg\n",&tilt[i],&angle[i],&udepth[i]) != EOF) 
i++; 
else notEOF=0; 
} 
if (LOWER_SENSOR==ON) 
{ 
j=1; /* read lower tilt data */ 
notEOF=1; 


while (notEOF) 


if (fscanf (fp2, "tg\n",&ldepth[j]) != EOF) j++; 
else notEOF=0; 

} /* while */ 

if (i<=j) LastDelay=i-1; 

else LastDelay=j-1; 


else LastDelay=i-1; 


/*** Set Delay Clock Times in Seconds ***/ 


delayClock [0] =0.0; 
for (i=1; i<=LastDelay; i++) delayClock[i] =(float)i*60.0; /* Delay 


Execution Times */ 


[/exxexxeeeeeeTHIS 18 the assumed array geometry: LINEAR*******eeex / 
for (j=1; j<=CHANNELS; j++) 


for (i=1; i<=LastDelay; i++) 
xx [J] [1] =DELTA_R* (float) (j-1) *cos (angle [i] /RADIAN) * 
sin(tilt [i] /RADIAN) ; 


yy [j] [1] =DELTA_R* (float) (j-1) *sin (angle [i] /RADIAN) * 
sin(tilt [i] /RADIAN) ; 
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zz[j] [1] =DELTA_R* (float) (j-1)*cos (tilt [1] /RADIAN) + 
OFFSET*cos (tilt [i] /RADIAN) +udepth [i] ; 
) aor, 
| / torre, 


[RR ek ieee kek ikke tee ae tote ate tee koe ate tee tee tek tek tek kk kk tk tek kk a I | 


theta=(360.0-LOOK_DIRECTION) /RADIAN; /* coordinate rotation */ 
for (j=l; j<=CHANNELS; j++) 


for (i=1; i<=LastDelay; i++) /* points x into signal */ 


x[j] (1) =xx[ 3] [1] *cos (theta) -yy[j] [1] *sin (theta) ; 
y (3) (i) =xx[j] [1] *sin (theta) +yy[j] [1] *cos (theta) ; 
Zaza ler 
} fF teres, 
\ o/s for4t/ 


/* Memory Deallocation */ 
1f (LOWER_SENSOR==ON) free ((char*) 1ldepth) ; 
free ((char*) tilt) ; 


free ((char*) angle) ; 

free ( (char*) udepth) ; 
free_matrix(xx, CHANNELS) ; 
free_matrix(yy, CHANNELS) ; 

free matrix(zz, CHANNELS) ; 

fclose (fpl) ; 

if (LOWER_SENSOR==ON) fclose(fp2) ; 
ng ehebbaeii( (i) 


[RRR KKK KK KR RK KR KR KR KK RK / 


/***** END processTilt *****/ 
[| 8 i ete ie te ke tke te te toe kee ite kk a / 


[xxeee 

* FUNCTION: processModes () 

* 

* This function handles the normal mode data. It calculates hydrophone 
* weights and group speed. The user must insure that the depth vector 
* and eigenfunction vector are of equal length. 

* 

* Arguments: z[] {) weight [] [] 
Dtrc 

* 

* Return value: 0 

* 

* Functions called: vector () ExitOnError () 
* INTERPOLATE () dydx () 

* 

* Definitions called: ANSI UNIX 

* PI SSP_LENGTH 

* EIGVAL_ LENGTH ORDER 

* CHANNELS F CARRIER 

* 

* Global variables called: Minute 

* 

* Significant memory allocation: depth [] Zn [] 

* OnoOfE [] w[] 

* Kr [] dwdK 

* 

kk ke / 


#if defined ( ANSI ) 
int processModes (float **z, float **weight, float *ptrC) 
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#elif defined ( UNIX ) 
processModes (z,weight, ptrC) 
float **z, **weight, *ptrcC; 
#Hendif 


int i, j, k, ptsEigVal, set, notEOF, deadPhones, weightNotAssigned; 

int items, moreData=ON; 

static int ptsDepth, ptsEigFun; 

Eigat “WwW, *Kige *“dwdK, *Workj;»eerr, df, freq; 

static float depth (SSP_LENGTH+1] BAY {[SSP_LENGTH+1] (EIGVAL LENGTH+1], 
OnOff [CHANNELS+1] ; 

char key, fileName [80] ; 

FILE *fpl, *fp2; 


if (firstBuffer==1) 
w= (float*) vector (EIGVAL_LENGTH) ; 
Kr= (float*) vector (EIGVAL_ LENGTH) ; 
dwdK= (float*) vector (E IGVAL_LENGTH) ; 
Work= (float*) vector (EIGVAL_LENGTH) ; 
printf ("Enter file name for normal mode data (eigenfunction): "); 


i1f((scanf("%$s", fileName) ) ==EOF) 
ExitOnError("Fatal error in scanf()"); 


Princt ("9 \n vn); 
if ((fpl=fopen(fileName, "rb")) == NULL) 
ExitOnError ("Error opening NORMAL MODE data file 
(eigenfunction) .") ; 


printf ("Enter file name for normal mode data (eigenvalues): "); 


if ((scanf("%s", fileName) ) ==EOF) 
ExitOnError("Fatal error in scanf()"); 


petntr ("\n\n") ; 


if ((fp2=fopen(fileName, "rt")) == NULL) 
ExitOnError ("Error opening NORMAL MODE data file 
(eigenvalues) .") ; 
i=1; /* read normal mode eigen values */ 
NOLcEOr=1 ; 


while (notEOF) 


if (fscanf(fp2, "sg sg \n", &w[i], &Kr[i]) != EOF) i++; 
else notEOF=0; 


ptsEigVal=i-1; 


/* read normal mode data Zn */ 
set= (ORDER+1) /2; 
j=1; 
df=2*PI* (float) F_SAMPLE/ (float) FFT LENGTH; 
while (moreData) 


items=fread (Work, sizeof (float) ,ptsEigVal+1,fp1) ; 


if (items==ptsEigVal+1) /* continue */; 
else if (ferror(fpl1) !=0) 


Wal 


ExitOnError("Error encountered while reading Eigenfunction 
File"); 
else if (feof(fpl1) !=0) moreData=OFF; /* EOF */ 
else ExitOnError ("Unknown error handling Eigenfunction data") ; 


depth [j] =Work [0] ; 

/* perform "quick" 2nd order polynomial interpolation between 
frequencies */ 

freq=df; 

k=1; 

i=1; 

set=0; 

while (freq<=w[ptsEigVal] ) 


if (freq>=w[1]) && ((float)i/((float)k+2) < (w[k+1]-wlk])/df || 
k+2 == ptsEigVal) ) 


fastpoly2 (&w[k-set] , &Work [k-set] , - (ORDER+1) ,freq, &Zn[j] [1] , &err) ; 
i++; 
fregq+=af; 
} /* if */ 
else if (w[k+1]<freq) k++; 
else freqt=df; 
} /* while frequency in range */ 
j++; 


/* while(moreData) */ 


fmin=w[1]) /(2*PI) ; 

fmax=w [ptsEigVal] /(2*PI) ; 

printf ("\nProcessing Bandwidth is from %g to %g Hz.\n\n", fmin, fmax) ; 
ptsEigFun=i-1; 

ptsDepth=j-1; 

for (i=1l;i<=ptsDepth;i++) depth[i] =depth[i]+CTD_OFFSET; 

for (i=1;i<=CHANNELS;i++) OnOff [i] =1.0; 

printf ("Do you want to turn off any hydrophones? "); 


if ((scanf ("%s", &key) ) ==EOF) 
ExitOnError ("Fatal error in scanf()") ; 


if (key==’y’ || key==’Y’) 
printf ("\nHow many hydrophones must be secured? ") ; 


if ((scanf("%i", &deadPhones) ) ==EOF) 
ExitOnError ("Fatal error in scanf()") ; 


for (i=1;1<=deadPhones ;i++) 


{ 


printf ("\nEnter hydrophone number to secure: "); 


if ((scanf ("%i", &j]) ) ==EOF) 
ExitOnError("Fatal error in scanf()") ; 


if (j>CHANNELS || j<1) 


ExitOnError ("Bad hydrophone identification") ; 
Guere lai -o- 0. 


TZ 





i Lor. * / 
comet #7, 


Ha aR aa 


for (i=1;i<=LastDelay;i++) 


if (depth [ptsDepth] <z [CHANNELS] [1] ) 


printf ("Max eigenfunction depth is: %f£\n",depth[ptsDepth] ) ; 
printf("at depth[] index of: %1\n",ptsDepth) ; 
printf ("Max depth of phone number %1 is: %f\n\n", 
CHANNELS, z [CHANNELS] [i] ) ; 
ExitOnError("Fatal data set error"); 


ere if */ 


jay * tor */ 


for (i=1; i<=CHANNELS; i++) 


He (OnOLt [i] ) 


for (k=1;k<=ptsEigFun;k++) 


weightNotAssigned=1; 
j=1; 
while (j<=ptsDepth && weightNotAssigned) 


if (z[i] [Minute]<0.0 || depth[j]<0.0) 


Breinee (<a \n" 21) ; 

Peimee (y—s1 \n", a)’; 

printf ("Minute=%i\n",Minute) ; 

pmimcn('2Z (1) [(Manutel is: t£\n", z([ij] (Minute]); 

printf ("depth[j] is: *f\n",depth[j]) ; 

printf ("Depth less than zero encountered in 


processModes.") ; 


orientation") ; 


z[(ij] [Minute], 


Beinet (Usn\n") ; 
ExitOnError ("Check input depths for coordinate 


has it  */ 


if (depth[j] <z[i] [Minute] && depth[j+1] >z [i] [Minute] ) 


set= (ORDER+1) /2; 
INTERPOLATE (&depth[j-set], &Zn[j-set] [(k], ORDER+1, 


&weight [i] [k], &err) ; 
weightNotAssigned=OFF ; 
est </ 
if (depth [j] ==z [i] [Minute] ) 
weight [i] [k] =Zn[j] [k]; 
weightNotAssigned=OFF; 
ease * / 


j++; 


/* while */ 


JX LOD © */, 
ee Ono et a / 
Ve / amcor i) */ 


We: 


/* set weights to zero if phone turned off */ 
for (i=1; 1i<=CHANNELS; i++) 


if (OnOff) /* do nothing */; 
else 
for (k=1; k<=ptsEigFun; k++) weight [1] [k] =OnOff [i] *weight [i] [k] ; 
\ </ * Sher 


if (Minute==1) 


dydx (Kr,w, dwdK, ptsEigVal) ; 
for (i=1;i<=ptsEigVal;i++) 


if (w[i]<2.0*PI*F CARRIER && w[i+1]>2.0*PI*F CARRIER) 


set= (ORDER+1) /2; 
INTERPOLATE (&w[i-set] , &dwdK[i-set] , ORDER+1, 
2.0*PI*F CARRIER, ptrc, &err) ; 
} /* if */ 
} /* for */ 


free ((char*) w) ; 
free ((char*) Kr) ; 
free ( (char*) dwdK) ; 
} /* if */ 
return( 0O ); 
KkKkkkkkkkkkkkKKKKKKkKKKKKkKKKKEK 
} / / 
[/***** END processModes *****/ 
[wR RR RR ik tt kK tk tk kik eK / 
[wwwx 


* FUNCTION: dydx() 

* 

* This function estimates derivatives. 

* 

* Arguments: x [] y{] 
* ddx [] points 
* 

* Return value: 0 

* 

* Functions called: ExitOnError () 

* 

* Definitions called: ANSI UNIX 
* STEP 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

* 

wk / 


#if defined ( ANSI ) 

int dydx(float *x, float *y, float *ddx, int points) 
#Helif defined ( UNIX ) 

dydx (x,y, ddx,points) 

float *x, *y, *ddx; 

int points; 

#Hendif 


int nn; 


for (n=1;n<=points;n++) 
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if ((n>=STEP) && (n<=points-STEP)) /*center*/ 
ddx [n}] = (y[n+STEP] -y [n-STEP] ) / (x [n+STEP] -x[n-STEP] ) ; 


else if (n<STEP) /*beginning*/ 
ddx [n] = (y [n+STEP] -y [1] ) / (x[n+STEP] -x[1]) ; 


else if (n>points-STEP) /*end*/ 
ddx [n] =(y [points] -y[n-STEP] ) / (x [points] -x[n-STEP]); 


else 
ExitOnError ("Index error in dydx"); /* sanity check */ 

je for */ 

meturn( 0 ); 
} [| & % ek a KK KK ee RK / 

[xxxwe END dydx tke kee / 

[| #8 ee ik ik kk tk kK tk tke / 
kkekkxk 
* FUNCTION: fastpoly2 () 
* 
* This function performs second order polynomial interpolation for 
interpolating 
* between modal frequency. 


* 

* Arguments: xa [] ya [] 
* n x 

. y dy 
* 

* Return value: 0 

* 

* Functions called: vector () ExitOnError () 
* 

* Definitions called: ANSI UNIX 
* 

* Global variables called: none 

* 

* Significant memory allocation: d[] eld 
* 

etek tee / 


#if defined ( ANSI ) 

int fastpoly2( float *xa, float *ya, int n, float x, float *y, float *dy) 
#elif defined ( UNIX ) 

fastpoly2 (xa,ya,n,x,y, dy) 

mGat *xa, “ya, X, *y, *day; 

anc en; 

#endif 


mc 1; 
float att, dx: 
if (abs(n) != 4) 
ExitOnError ("USER error, fastpoly3 interpolater can only be ORDER 
2) ee 
dx = xa[1]-xa[0]; 
if (n == 4) /* FIND CLOSEST POINT */ 
if (x <= xa[3]+dx/2 && x > xa[3]-dx/2) i=2; 


else if (x <= xa[3]-dx/2 && x > xa[2]-dx/2) i=1; 
else if (x <= xa[2]-dx/2 && x > xa[1]-dx/2) i=0; 


iS 


else 
ExitOnError ("Error in routine FASTPOLY2") ; 


else i=1; 

dif=x-xa [i]; 

“Vaya vale 

*y-= Gif*(3.*ya[i] - 4.*ya[i+1] + ya[i+2]) /(2.*dx) ; 
*y+= dif*dif*(yali] - 2.*ya[i+1] + yali+2]) /(2.*dx*dx) ; 


*dy=0; /* error term to be installed later */ 


return (0) ; 


} [| % Hi ie te te He ek te te te tee eke kk tek / 


[exe END fastpoly2 elielelelell | 


[ei ie te ee ee eee tee te kee eee Heke / 


[**ewe 

* FUNCTION: fastpoly3 () 

* 

* This function performs Third Order polynomial interpolation. 
* 

* Arguments: xa [] yal] 

* n x 

i y dy3 

* 

* Return value: 0 

¥ 

* Functions called: vector () ExitOnError () 
* 

* Definitions called: ANSI UNIX 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

¥e 

te te ke de / 


#if defined ( ANSI ) 

int fastpoly3( float *xa, float *ya, int n, float x, float *y, float *dy3) 
#elif defined ( UNIX ) 

fastpoly3 (xa,ya,n,x,y,ay3) 

ficat *xa,; *ya, x; *V7e2cy2; 

int ne 

#tendif 


pat 1: 
flioat Gif, dx yo, 


if (abs(n) != 4) 
ExitOnError ("USER error, fastpoly3 interpolater can only be ORDER 
3 
Sa 


ax = xa[iti1] -xal[i] ; 

@dif =x - xa la-i); 
yo= 4*ya[1] - 6*yal2] + 4*ya[3] - ya[4]; 
*dy3 = -yall] + 3.*yal2] - 3.*ya[3] + ya[4]; 


*y = yo + (dif/dx*( -13./3.*ya[1] + 9.5*ya[2] -7.*ya[3] + 11./6.*yal4]+ 
(dii/(2.*Gx)* ( a yvaii) 8.*ya[2] + 7.*ya[3] - 2.*yal4] + 


vc 


(ites. *dx)*( *dy3 )))))); 


return (0) ; 


} I heainatietieiatietaieatiahatetehetehetaheteheteheiehel 
/*x*** END fastpoly3 *****/ 
[#4 i Hee ie ee ie ee tee tee te tee te ie kee ie tee / 
[ewww 


* 


% 


FUNCTION: polint () 


* This function performs polynomial interpolation. 

* 

* Arguments: xa [] ya[] 
* n x 

id y dy 

* 

* Return value: 0 

* 

* Functions called: vector () ExitOnError () 
* 

* Definitions called: ANSI UNIX 
* 

* Global variables called: none 

* 

* Significant memory allocation: da{[] et] 
* 

wk ak / 


#if defined ( ANSI ) 

int polint( float *xa, float *ya, int n, float x, float *y, float *dy) 
#elif defined ( UNIX ) 

polint (xa,ya,n,x,y,dy) 

float *xa, *ya, x, *y, *dy; 

int n; 

#Hendif 


int 1, m, ns=1; 
float den, dif, dift, ho, hp, w; 
Eloat *%c, *d; 


n=abs (n) ; 
dif=fabs (x-xa[1]) ; 


c=(float*) vector (n) ; 
d=(float*) vector (n) ; 


for (i=l; i<=n; i++) 
if ((dift=fabs(x-xa[i])) < dif) 


ns=1; 
dif=dift; 
} /* if */ 
c[i] =yal[il ; 
d[i} =ya [i] , 
fe) om ~/ 
*v=ya [ns--] Ps 
for (m=1; m<n; m++) 


for (i=l; i<=n-m; i++) 


{ 


ay 


ho=xa [1] -x; 
hp=xa [1+m] -x; 
w=c [1+1] -d[i] ; 
if ((den=ho-hp) ==0.0) 
ExitOnError ("Error in routine POLINT") ; 
den=w/den; 
d[i] =hp*den; 
c [i] =ho*den; 


} stor S, 
aes he a < (n-m) ? c{ns+1) : d[ns--])); 
* or * 


free ((char*)d) ; 
free((char*)c); 
return( O ); 


} [#7 ee ee ie te de te ie ee ee ke tee tee ee / 


[xe ka END polint wie tee / 


[RRR Ka RK aK kk kk ke / 


[xxkKe 


* 


% 


FUNCTION: rataint () 


* This function performs rational function interpolation. 
* 

* Arguments: xa [] ya [J 
a n x 

ui y dy 
* 

* Return value: 0 

* 

* Functions called: vector () ExitOnError () 
* 

* Definitions called: ANSI UNIX 
* TINY 

* 

* Global variables called: none 

* 

* Significant memory allocation: d[] c[] 
* 

ee te kee / 


#if defined ( ANSI ) 

int ratint( float *xa, float *ya, int n, float x, float *y, float dy) 
#elif defined ( UNIX ) 

ratint (xa, ya,n,x,y, dy) 

float *xa, *ya, x, *y, *dy; 

int n; 

#Hendif 


{ 


Int Mm, 2; NS=l> 
float Wer Gs; hh, h, dd, ey *<c; 


n=abs (n) ; 

c= (float*) vector (n) ; 
d= (float*) vector (n) ; 
hh=fabs (x-xa[1]) ; 
for (i=1; i<=n; 1++) 


h=fabs (x-xa[i]); 


if (h==0.0) 
{ 
*ysya [i]; 
*dy=0 .0 7 


free ((char*) d) ; 


is 


free ((char*)c) ; 
return ( 0 ); 

} /* if */ 

else if (h<hh) 


ns=i; 
hh=h; 

} /* else if */ 

c(i] =yal[i] ; 

d[i] =ya [1] +TINY; 
ieee tor * / 
*ysyal(ns--]; 
for (m=1;m<n;m++) 


for (i=1; i<=n-m;1i++) 


w=C [141] -d[1] ; 
h=xa [i+m] -x; 
t=(xa[i] -x)*d[i]} /h; 
dd=t-c[i+il] ; 
if (dd==0.0) 
ExitOnError ("Error in routine RATINT") ; 


dd=w/dd; 
dad [i] =c[1i+1] *dd; 
c [i] =t*dd; 
faye for */ 
’ ie Ss < (n-m) ? c[ns+1]) : d[ns--])); 
or * 


free ((char*)d) ; 
free ( (char*)c) ; 
return( 0 ); 

} [RRKKKKKKKK KK KKK KKK KKK / 
/**xx** END ratint ****/ 
[RRR RK KR KR KR KK RR RK KK / 

[eK 


* FUNCTION: window () 

* 

* This function applies a Blackman window to a vector. 
* 

* Arguments: data [] N 
* 

* Return value: 0 

* 

* Functions called: none 

* 

* Definitions called: ANSI UNIX 
k PI 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

* 

wt ke ee / 


#if defined ( ANSI ) 

int window(float *data, int N); 
#Helif defined ( UNIX ) 

window( data, N ) 

float *data; 

Pit oN; 

tendif 


Us, 


int “a 


for (n=0; n<N; n++) 


data [n+1] =data [n+1] * (0.42+0.5*cos (2.0*PI* (float) (n-N/2) / (float) (N-1) ) 
+0.08*cos (4.0*PI* (float) (n-N/2) / (float) (N-1))); 


} /* foes, 
return ( 0 ); 
} [8 8 i i te te te te te keke ke ke tek tek tok / 


[eK KK END window tk tee / 
[| 8 iH i te te te de tee tek tek te toe kek / 


[® wR Re 

* FUNCTION: realft () 

* 

* This function calculates FFT’s 

* 

* Arguments: data [] 
* isign 
* 

* Return value: 0 

rd 

* Functions called: fourl () 
% 

* Definitions called: ANSI 
* 

* Global variables called: none 
* 

* Significant memory allocation: none 
% 

wk ke tet / 


#if defined ( ANSI ) 

int realft (float *data, int n, int isign) 
#elif defined ( UNIX ) 

realft (data, n, isign) 

float *data; 

int ne Sloan: 

#endif 


int i,> Pl; 1251356147 en2ps: 
float cl=0.5, €2, hixr, bli, H2r,- hea 
double wr, Wi, wpr, wpi, wtemp, theta; 


theta=3 .141592653589793/ (double)n; 
lf (isign==1) 


c2 = -0.5; 

fourl (data, n, 1); 
} | Coe a * / 
else 

c2=0.5; 

theta = -theta; 


} /* else */ 

wtemp=sin (0.5*theta) ; 
r= -2.0*wtemp*wtemp; 

wpi=sin (theta) ; 

wr=1.0+wpr; 

wi=wpi; 

n2p3=27e35), 
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UNIX 


for (1=2; i<=n/2; i++) 

14=1+ (13=n2p3- (12=1+ (11=i1+i-1))); 
hir=cl1* (data [i1]+data[i3]) ; 
hli=cl1* (data [i2] -data[i4] ) ; 
h2r = -c2*(data[i2]+data[i4]); 
h2i=c2* (data[il] -data[i3]); 
data [il] =hlr+wr*h2r-wi*h2i; 
data [i2] =hli+wr*h2i+wi*h2r; 
data [13] =hlr-wr*h2r+wi*h2i; 
data[i4] = -hli+wr*h2i¢+wi*h2r; 
wr= (wtemp=wr) *wpr-wi*wpi+wr; 
wi=wi*wpr+wtemp*wpitwi; 

es *® for */ 

if (isign==1) 


data [1] =(hlr=data[1] )+data[2] ; 
data [2] =hlr-data[2] ; 
yt it */ 


else 


data [1] =c1* ((hlr=data[1])+data[2]) ; 
data [2] =c1* (hir-data[2]); 
fourl (data,n,-1); 
} /* else */ 
fecurm ( 0 ); 
} [8 RR Rae te kek kk tok dei ieee / 
[xxxw* end realft ielielieieliell | 
| cheiichehehioheiehelohehelshetehelehelehehehel | 
[xxx 


* FUNCTION: fourl () 

* 

* This function calculates FFT’s 

* 

* Arguments: data [] nn 
* isign 

* 

* Return value: 0 

* 

* Functions called: none 

* 

* Definitions called: ANSI UNIX 
* SWAP 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

* 

wk ke / 


#if defined ( ANSI ) 
mitetourl(float “data, int nn, int isign) 
#Helif defined ( UNIX ) 
fourl (data, nn, isign) 
float *data; 
tate, LSign,; 
#Hendif 
int n, mmax, Mm, j, istepyw i; 


double wtemp, wr, wpr, wpi, wi, theta; 
float tempr, tempi; 


idl 


N=nn <<; 


j=l 


e 
# 


for (fal sic 14-2) 


nl! 


LE eee) 


SWAP (data[j] ,data[i]) ; 
SWAP (data [j+1]) ,data[i+1]) ; 
) Ve ecrss, 
M=n >oe1 
while (m >= 2 && j > m) 


jJ -= mM; 
Moose 
} /* while */ 
5 += Mm; 
* for */ 


mmax=2 ; 


whi 


le (n > mmax) 


istep=2*mmax; 

theta=6 .28318530717959/ (isign*mmax) ; 
wtemp=sin (0.5*theta) ; 

wpr = -2.0*wtemp*wtemp; 

wpi=sin (theta) ; 

wr=1.0; 

wiz=0.0; 

for (m=1; m<mmax; m+=2) 


for (i=m; i<=n; i+=istep) 


j=1+mmax; 
tempr=wr*data [j] -wi*data [j+1] ; 
tempi=wr*data [j+1] +wi*data[j] ; 
data [j] =data [i] -tempr; 
data [j+1] =data[i+1] -tempi; 
data[i] += tempr; 
data[i+1] += tempi; 
oey ha asong 2) 
wr= (wtemp=wr) *wpr-wi*wpit+wr; 
wi=wi*wpr+wtemp*wpitwi ; 
} /* for 4% 
mmax=listep; 


} /* while */ 


ree 
kkk 


man ( 0 ); 
kkk kkk kkk kK / 


[RRR end fourl kk kk / 


V hohokhokhotehahehehoheihehehahehehehelel | 
[ke KKe 
* FUNCTION: vector () 
* 
* This function allocates memory for UNIT OFFSET vectors. 
* 
* Arguments: length 
* 
* Return value: wis 
* 
* Functions called: ExitOnError () 
* 
* Definitions called: ANSI UNIX 
* 
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* Global variables called: none 
* 

* Significant memory allocation: v [] 
* 

KkKKKK 


#if defined ( ANSI ) 
float *vector(int length) 
#elif defined ( UNIX ) 
vector (length) 

int length; 

#endif 


float *v; 


if ((v=(float*) malloc ((length+1) *sizeof (float) ) ) ==NULL) 
ExitOnError ("Memory allocation failure in vector().") ; 
#if defined ( ANSI ) 
return v; 
#elif defined ( UNIX ) 
return (long int)v; 
#Hendif 


He He de ee ke ae de Ye ke oe ae ae te ke ae ae ae te ke ie ae / 


[/*x*x*** END vector ete He kee / 
[| % He He He He ee ae te te ke He ae eke ke He He / 


[exw 

4 FUNCTION: matrix () 

* This function allocates memory for UNIT OFFSET 2-D arrays. 
* 

* Arguments: row col 
* 

* Return value: * eM 

* 

* Functions called: ExitOnError () 

* 

* Definitions called: ANSI UNIX 
* 

* Global variables called: none 

* 

* Significant memory allocation: mf) [J 

-... 


#if defined ( ANSI ) 

float **matrix(int row, int col) 
#Helif defined ( UNIX ) 

matrix (row,col) 

iat row, col; 

#Hendif 


at 1: 
float **m; 


if ((m=(float**)malloc( (unsigned) (row+1) *sizeof (float*) ) ) ==NULL) 
ExitOnError ("Allocation failure 1 in matrix()"); 
for (i=l; i<=row; i++) 


if ((m[1i] =(float*)malloc( (unsigned) (col+1) *sizeof (float) ) ) ==NULL) 
ExitOnError("Allocation failure 2 in matrix()"); 
/* for */ 
#if defined ( ANSI ) 
return mM; 
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#Helif defined ( UNIX ) 
return (long int)m; 
#endif 
| ch chshahehehehiehehehehehehehehehehehehehel | 
[RRKKK END matrix wk ake / 
hee te He te te ke te ke te te ke tee ie kok ke te tet te / 
te ke 4 


* FUNCTION: free_matrix() 


* 

* This function deallocates memory from UNIT OFFSET 2-D arrays. 
4 Arguments: m{j [J] row 
* Return value: 0 

* 

* Functions called: none 

* 

‘ Definitions called: ANSI UNIX 

* Global variables called: none 

* 

* Significant memory allocation: none 

a. 


#if defined ( ANSI ) 

void free matrix(float **m, int row) 
#elif defined ( UNIX ) 
free_matrix(m, row) 

Eloate =~m); 

int row; 

#endif 


inte; 


for (i=row; i>=1; i--) 
free ((char*)m[i]J) ; 
free ((char*)m) ; 
return( 0 ); 
I ishahahchehstehehshelehehehehehehehehehehehehehehehekel f 


/***** END free matrix *****/ 
[8 kk kk ite te te te te ee ek ke kk Ke / 


[eR RK 

* FUNCTION: ExitOnError () 

* 

* This function performs an abnormal process termination. 
* 

* Arguments: error txt [] 

* 

* Return value: none 

‘ Functions called: none 

* 

* Definitions called: ANSI UNIX 
‘ Global variables called: none 

: Significant memory allocation: none 

ees, 


#if defined ( ANSI ) 
void ExitOnError(char error_txt []) 
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#elif defined ( UNIX ) 
ExitOnError (error txt) 
char error txt []; 


#endif 
Eprantt (stderr, "Program run-time error ...\n"); 
fprintf(stderr, "ts\n",error_ txt) ; 
fprintf (stderr,"...now exiting to system...\n") ; 


fclose (fpInSound) ; 

if (SIGNAL==ON) fclose (fpOutSound) ; 

if (ERROR_ESTIMATE==ON) fclose(fpMSE) ; 
exit (0); 


2h ee te te te tee a ke ae ie te a ke ae ae eke eke ee ee / 


/***** END ExitOnError *****/ 
[#8 RRR ek ke a kk kk kk ee / 


B. MODAL DECOMPOSITION PROGRAM 

This program performs modal decomposition for a specified 
bandwidth and frequency increment. Inputs to the program 
include sound speed and density profiles. An example input 
file can be found at section D of this Appendix. This program 
also uses a complex eigen solver called as an external 
subroutine. Output from this program is a direct access 
binary file containing the eigen solutions evaluated at the 


specified frequencies. 


* Last Updated by G.A. Omans II 
* enesO uuly 1992 


* 
kaekkkkkkkkkkkekKkKkKkKkKKKKKEKKEKKKKKKKKKEKKKKK KK KKK KKKKKKKKKKKKKKKKKKKKKKKKKKKEK 
kkkk« 
* MODAL DECOMPOSITION PROGRAM 
* 
Program to calculate eigenmodes and eigen values given the sound 

& density profile 

Output is to a DIRECT ACCESS Binary file "modes. dabin" 

Program used to extract Neccesary Data for Beam Forming is 

"CWextract.f" 


t + + F 


* or "BBextract.f". CWextract creates input to the Continuous-Wave 
* beamformer. BBextract creates input files for use in the Broadband 
* beamformer. 


kkekkkkkkkkkkkkKKekKkkKKKKKKKKKKKKKKKKRKKKK KKK KKK KKK KKK RK KK KR KK aK a RK KK kK 


kkk kk 

* FILES: : 

* afreqmode.in - Input file containing sound & density profiles, 
* depth increment dz, hard bottom depth Hdepth, 

* center frequency fc (Not used in computations) . 


3) 5 


* £ * &€ & Fee EF EE F FE Ee FE FE EE HK KE KK HE EE HF FE HE FH OF 


Profiles need not be complete, they will be inter- 
polated to correct number of points. However, the 
profile depth spacing must be greater than dz at all 
points. 


modes.dabin - Direct access, binary output file. This file can be 


accessed by "mextract.f" to extract modal information 
for beamforming, or for inspecting mode shapes & 
eigenvalues. 


mode .sys - This file is used as a check on how well the eigen 


solver "eigrf" and the eigen subroutine itself are 
doing. A fair amount of information is output to 
this file, including sound & density profiles, the 
eigen matrix a, the number of trapped modes & their 
eigenvalues; to name a few. If VERBOSE parameter 
option is selected, then even more information is 
displayed. 


standard output - Information sent to standard output is very 


imited. 


The user is provided with just enough information to 
indicate where the program is currently running, and 
any Significant problems if encountered. If you 
have selected VERBOSE prior to compiling the program, 
then the amount of information will overwelmingly 
increase, and will be impossible to keep track of in 
the interactive mode. In this case the program should 
be run in the background and standard output should 
be redirected to a file. In any case, the program 
has a significant run time do to the recursive eigen 
solver routine, and should be run in the background 
most of the time. 


kkekkkkkkkkkkkkkkkkkkkkkkkkkkekkkkkkkkkkkkekekekkkekkkkekkekk keke KKK KKKKKRKKRKKKKK 
kkk kk 

*QUIRKS : 
Because of the recursive operations performed by IMSL subroutine 


t+ £ + & & F 


elgrf, 


a large amount of memory is required. To accomadate the 


memory needs of the subroutine, in the main program the arrays Z 
and kc have been demensioned larger than required by eigrf itself. 
This extra cusion of memory significantly improves the solutions 


from the eigen solver and increases the number of trapped modes found. 
KkKkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk kkk kkk kkk kkk kkkk kkk kkk 


kkekkk 


* * &£ + € € FF Fe SF FF FF SF KE OF 


a - 
buee - 


Cc 
depth = 


ebe - 
dz z 
efun = 
eigen - 
EX = 
- - 
EG = 


fdate 
filename 


(array) eigenvalue coefficients from difference eqn. 
(small dz approx.) 
(array) profile buffer for input into READP subroutine 
(array) sound speed, interpolated from input file 
(par) used to define the soft depth; This parameter 
does not effect calculations, however, it is used for 
testing the reliability of trapped eigen solutions at 
the silt. 
(par) frequency increment used in evaluating band 
(var) depth step read from input file 
(array) eigen functions 
SUBROUTINE to calculate eigenfunctions and eigen values 
(log) switch used for testing the existance of a file 
(var) frequency (Hz) 
(var) center frequency; serves no purpose except as identifier 
- (char) used for printing start and end times 
- (char) holds the name of a file for error messages 
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* fmax - (par) maximum frequency in band to evaluate 

* fmin - (par) minimum frequency in band to evaluate 

* Hdepth - (var) depth of channel + the sediment depth 

* le. depth of hard bottom 

* ifreq - (var) frequency counter 

= icount - (var) number of terms read from profile for input into buff 
* isize - (parameter) size of temporary array used for interpolating 
profile 

* NOTE: isize must be > IFIX(1.5 + ZMAX profile/DZ) 

* irec - (var) record number for direct access file 

* itrap - (var) passed by subroutine eigen indicates number of trap modes 
found 

* jump - (par) jump is the maximum number of "Trapped" Modes DESIRED 
= jump is also used to jump to certain records in the direct 
access 

* fille. 

=k - (array) eigenvalues or the square of Kn for each mode 

~ Wi - (par) number of modes to be calculated; for best results, 

* this should be isize-1l 

a nt - (var) total number of frequencies to calculate nf = 
1+ (fmax-fmin) /df 

* nmodes - total number of eigen solutions for a particular frequency 
written 

ie to direct access output file. (minimum between itrap & jump) 
= nz - (var) number of eigenvalues, profile points, & eigenvalue 
length 

- nZp1 - (var) number of points in svp (includes surface point, which 
is 

* not required for eigenvalue since the value 1s zero at the 
surface) 

* pi - (var) pi itself in double precision 

Salem k - (var) record length of direct access file (important to 
mextract.f) 

* rho - (array) density array,contained in input file mode.in 

* temp - (array) a working array used in subroutines readp and 


eigen; in 
* @igen temp is called mstor 


* ow - (var) radial frequency (rad/s) 

* VERBOSE - (log) Gives User full Status Report on All Eigen Solutions 
* and Profiles; Good For Trouble Shooting But Maybe 
Overwhelming 

cewek kekckekeKkKeKkeK KK KKK KKK KE KK KK KKK KKKKKKKKKKKK KK KKK KK KKKKK KK KKK KK KKK KKK KKKRKKKEK 
kkekek«x 

* NOTE: The input profiles are double precision. 

* Profiles are changed to double precision in READP subroutine. 
* Profiles in the input file need not be in double precision. 
eke mK KKK KKK KKKKEKKK KK KK KKK KKK KKK KKK KKK KKKKKKKKKKKKKKK KK KK KKK 
kekkx 

c 

Cc SYSTEM PARAMETERS 

Cc 


real*8 df,depth 
logical ON,OFF, VERBOSE 
parameter (ON=1, OFF=0) 


USER PARMAMETERS 


AAD 


parameter (VERBOSE=OFF) 

parameter (jump=20) 

parameter (f£min=208, fmax=240, df=1) 
parameter (isize=211,m=(isize-1)) 
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parameter (depth=276.) 


implicit real*8 (a-h,o-z) 

real*8 k(m),w,f,fc 

real*8 Hdepth, dz 

real*8 pi 

integer irec,nmodes,rcl 

character*24 fdate, filename 

logical EX 

DIMENSION c(isize) ,rho(isize) , temp (isize) 
DIMENSION buff (isize, 2) 


Cc 
@: EIGEN SUBROUTINE VARIABLES 
Cc 
DIMENSION ea(isize,isize) , work (2* (isize+2) *isize) 
DIMENSION kc (isize+40) ,z(isize+40,1isize+40) ,efun(isize,m) 
Cc 
Cc FILES 
C 
OPEN (UNIT=1, STATUS=’ OLD’ , FILE=’ afreqmode.in’ ) 
c 
Cc ese see eweeeeeeeeeeeeeeeeeeeeeeeeeereeeeeeeeeweewe ee @= @= 
c MORE open statements 
ee eS OE ee es 
e 
CG OPEN MODE.SYS 
e 
INQUIRE (FILE=’mode.sys’ , EXIST=EX) 
IF (EX) THEN 
OPEN (UNIT=13, FILE=’mode.sys’ , STATUS=’ OLD’ ) 
CLOSE (13, STATUS=’ DELETE’ ) 
ENDIF 
OPEN (UNIT=9, FILE=’mode.sys’, 
1 FORM=’ FORMATTED’ , STATUS=’ NEW’ , ERR=2003) 
ce 
c OPEN MODES .DABIN 
G 
nf = 1+((fmax-fmin) /df) 
IF (m.LE.isize.AND. (nf+2) .LE.isize) then 
rcl=isize*8 + 8*21 
ELSEIF ((nf+2).LE.m) then 
rcel=m*8 + 8*21 
ELSE 
rcl=(nf+2)*8 + 8*21 
ENDIF 
INQUIRE (FILE=’modes.dabin’ , EXIST=EX) 
IF (EX) THEN 
OPEN (UNIT=14, FILE=’modes.dabin’ , STATUS=’ OLD’ ) 
CLOSE (14, STATUS=' DELETE’ ) 
ENDIF 
OPEN (UNIT=4, FILE=’ modes .dabin’ , ACCESS=’ DIRECT’ , RECL=rcl, 
1 FORM=’ UNFORMATTED’ , STATUS=’ NEW’ , ERR=2005) 
& 
Cc DECLARE PI IN DOUBLE PRECISION 
c 
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gaa aaQ ree) qgqagqan aaAN PD lo 


QAAQ ANN 


QAaAN 


pi = 4.d0 * datan(1.d0) 


Give START TIME 


write(6,*) fdate() 
write (6,*)‘’\n’ 


WRITE (6,*) ‘'RECORD LENGTH: ‘’,rcl 


INITIALIZE BUFFER VARIABLE 


do 10 i=1,18s1ze 
do 20 j=1,isize 
bweet (7), 1) = 0.0 
20 continue 
10 continue 


Read in ocean/model parameters 


read(1,*) Hdepth 

write (6,*) ‘depth w/ sediment = ’,Hdepth 
read(1,*) dz 

read(1,*) fc 

write(6,*)'dz =',dz 

write(6,*) ‘frequency (Hz) = ',fc 


DETERMINE NUMBER OF POINTS IN PROFILES 
nzpl = 1.5 + Hdepth/dz 


aA = )elea| s)s Real 


Read in and interpolate sound velocity profile 
Read Sound Velocity Profile into Buffer Variable & Print 


i = 0 
do 30 WHILE (buff (1,1).NE.-1 .AND.1i1.LT.1i8s1zZe) 
1 = 141 
read(1,*) (buff(1,3) ,j=1,2) 
30 continue 
icount = i-1 
write(6,*) ‘SVP icount = ’,1icount 


Read and interpolate SVP 
CALL READP(dz,nzp1,isize,c,temp, buff) 
CLEAR BUFFER 
do i=1,isize 
bute ia) =O 
buff (1,2) =0 
enddo 


Write the interpolated SVP to the screen 


if (VERBOSE) then 
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OG:Q 


Q OY OOO Slee. Q OOo 


Q 
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do 60 j=1, (nzpl1+1) 
write (6,*) dz*(j-1), c(3) 
60 continue 
endif 


READ IN DENSITY PROFILE 


1=0 
DO 31 WHILE (buff (1,2) .NE.-1 .AND.1i.LT.isize) 
1 = itl 
read(1,*) (buff (i,j) ,j=1,2) 
31 continue 
icount = i-l 
write(6,*) ‘Density icount = ’,icount 


Read and interpolate DENSITY 
CALL READP(dz,nzp1,1isize, rho, temp, buff) 
Write the interpolated DENSITY to the screen 


if (VERBOSE) then 
do 61 j=1, (nzp1+1) 
write (6,*) dz*(j-1), rho(j) 
61 continue 
endif 


CLEAR WORKING ARRAY TEMP FOR USE AS MSTOR IN EIGEN SUBROUTINE 


do j=1,isize 
temp(j) = 0 
enddo 


CALCULATE EIGEN FUNCTIONS & EIGEN VALUES 
FOR ALL FREQUENCIES 


** nf = 1+((fmax-fmin) /df) ** CALCULATION PRIOR TO OPENING DATA FILE 
f = fmin 
write (6,*)'nE ="samm. di = der 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkekkkkkKkkkk KKK 


WRITE (4,REC=1) fc,fmin,fmax,df,nf,dz,m,nzp1,jump,rcl 
WRITE (4, REC=2) c. 
WRITE (4,REC=3) rho 


kaekkkkkkkkkkkkkkkkkkkkkkkkkkkkekkkkekekkekkkkkkkkkkkkkk kk kK 


irec=4 


do 100.1 freq = 91 nL 


CONVERT TO RADIAL FREQUENCY 
w=2* pi*f 
write (6, *) 2 \m\m¥ & % kk ae tee de ede dee tke de ke kde eke dete te detect te be tote tee tote teke ke! 


write(6é,*) ’f = ne (Hz) ifreq =',ifreq 


90 


write(6,*) ‘w ='’,w,’ (rad/s)’ 


INITIALIZE VARIABLES 


QAAN 


do i=l,nz-1 
k (i) =0 
ke (1) =0 
do j=1,nz-1 
i1f(j.LT.m) efun(i,3)=0 
z(i,j)=0 
enddo 
enddo 


CALL EIGEN SOLVER SUBROUTINE 


QAAN 


CALL eigen(c,rho,dz,nz,efun,m,k,kc,w,z,ea, 
& depth,work,temp,ifreq,isize, 
& itrap, VERBOSE) 


Cc kkk&kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


write (4,REC=irec) w,itrap,k 
ro KwKeKKKKKKKKKKKKKKKKKKKKKKKEKKEKKEKKKKKKKKKKKeKKKKKKRKKK 


if (jump.GT.itrap) then 
nmodes = itrap 

else 
nmodes = jump 

endif 


do 70 im=1, (nmodes) 
1rec=irec+l 


c KkKkkk&kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkekkk 


write (4,REC=irec) (efun(i,im) ,i=1,nz) 
ro. Ke KK KKK KKKKKEKKKKKKKKKKKEKKKKKKKKKKKRKKKKKK KKK KK KK 


70 continue 
irec=(jump+2) *ifreq + 4 


100 f=f+df 


go to 2500 
C eweeeeee ewe eweeweeeevee eewreeewerwrrFere rere wvweeewre wf=ere Freee wKee ws weweewrewvreweeee = = 
c close statements 
ee 
2003 filename=’mode.sys’ 
goto 2010 
2005 filename=’modes.dabin’ 
2010 write (6,3000) filename 
GOTO 2500 
2500 close (9) 
close (4) 
close (1) 
Cc 
Cc FORMAT STATEMENTS 
Cc 


eA. 


3000 format (’ ERROR OPENING FILE : ‘,A) 


END TIME 


write (6,*)’Program End Time: ‘',fdate() 


stop 
end 


kKkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkekkkkkkkKkkkkKkkKkkKK 
kk kkk 
KkKk&kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkKkkkkkkkkkkke 
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Cc 


SUBROUTINE EIGEN 


KRKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK 
kkk kk 


* 


+ 


+ ¢ €¢ + + € + + € FF FF FE FE OF FO OF OF OH OF 


INCLUDING DENSITY AND SOUND SPEED VARIATION 
NORMALIZED 
:eigenvalues,k and eigenvectors,z are complex 


compute m normal modes z(i,m) and wavenumbers k(1) 
for a given sound-speed profile c(i), 

a given frequency w, 

a given density profile 

and a given step size dz. 


nz*dz is the depth of the ocean 
be: pressure release surface. 
rigid bottom with soft sediment interface 
input: nz,m,dz,w,c 
output: w 
m,nz 
c,k,z (freq.,sound speed,eigenvalues, eigenvectors) 
(note: eigenfunctions,z have two less points than ¢; 
recall z(0) is O and z(nz)=z(nz-1). 
internal: a,work; one must set iw=15*n 


Kk&kKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK KKK KKK KKK KK KK KKK K 
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UPDATES /CORRECTIONS : 


1) UPDATED - OUTPUT ONLY TRAPPED MODES BASED ON 

K AT THE BOTTOM & ELIMINATES EVANESCENT 
SOLUTIONS. 15 JUN 92 

2) CORRECTION - EIGENMODE NORMALIZATION NOW REFLECTS 
INVERSE OF THE DENSITY VS. DENSITY MULTIPLICATION. 
18 JULY 92 

3) UPDATED - REORDER EIGEN SOLUTIONS BY MODE NUMBER. 
REORDERING IS ACCOMPLISHED BY SORTING BY EIGENVALUE 
FROM HIGHEST TO LOWEST. THIS WILL EFFECTIVELY ORDER 
ALL MODES FROM LOWEST TO HIGHEST. 24 JULY 92 

4) UPDATED/CORRECTED - REALIGNS TRAPPED MODES TO HAVE THE 
SAME SIGN OR PHASE. REFERENCE IS A NEGATIVE SIGN AT 
THE BOTTOM. 29 JULY 92 

5) UPDATED - CHECKS FOR BAD EIGENFUNCTION SOLUTIONS BY 


EVALUATING THE EXISTANCE OF MODAL ENERGY AT THE BOTTOM. 
30 JULY 92 


92 


* 


KeKKKKKKKEKKEKEKKEKKEKEKKEEKKEKKAKRKKEKKKKEKKKKKKKRKKKEKKKKaKKKKKKKKKKKKKKKKKKKKK KKK 


* 
* 
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calculating 


* 


* * +1 + + € € + + FF F EF EE € € FF KF OF OF OF KE OF KF OH OF OF 


* 
* 
* 


a 
all 
S 
eb 


d 
depth 
dterm 


dz 
dz2 


dz2inv 
dz2inv2 


efun 
eigrf 
ext 
itrap 


17Ob 


1k 
1size 
jlast 
k 

kb2 
KC 


k_ trap 
k_last 


m 


mode.sys 


mreal 
mstor 


nZ 


ize 


nzpl 
pi 


VERBOSE 


W 
work 
Zz 


finite difference matrix or "eigen matrix" 
- integer defining domain of eigen solutions 
sound speed profile for an individual station 
- value of the sound speed at the bottom of the water column 
sediment 
density profile for an individual station 
- depth of soft bottom, used to test for energy at bottom 
- a term used in the finite difference equation for 


the eigen matrix 
- depth increment of input 
- dz squared, used in finite difference EQ 
inverse of dz2 
- two times the inverse of dz2 
- eigenfunction matrix (Eigen Vectors) Represents Modes 
- IMSL matrix eigenvalue solver subroutine 
- Minimum value for sign testing of eigen function 
- counter used in conjunction with mstor to save only 
trapped modes (FINAL itrap IS TOTAL TRAPPED MODES) 
- =i1: calculate k 
2: calculate k & z 
3: calculate k,z & performance index 
- simular use as im 
- used as a modal incrementer 
- a demensioning term 
- integer variable used to check number of modes realigned 
real eigenvalue vector, k = Kn**2 
- eigenvalue at the bottom of water column sediment squared 
- complex eigenvalue matrix 
variable used to relay trapped Kn’s 
test variable used in reordering eigen solutions 
number of modes Evaluated from eigrf output 
- System file for subroutine, make sure you check "p index" 
- Number of Modes Realigned (CURRENTLY DOESN’T WORK) 
- array used to keep track of only modes which are trapped 
efun are then stored as a function of mstor 
- nzpl-1, number of depth points; also is eigen function 


- number of physical points (including free surface) 
- itself pi, in double precision 
- logical switch for print out control 


- frequency (2*pi*f) 


- internal work space: require for IMSL routine "eigrf" 


- eigenfunction matrix output from eigrf subroutine 
Kk kee aKa KaKKaKaKKKKKKKKKKKKKK KKK KK KaKK KK KKK KKK KKK RK RK KK KR KR KR IK KK KG a KK kK 


KH KKK KKH HKHKAKEKKAKKKKRKKKKKKKKKKKKKKKKKKRKKR KK KKKRKKKKKKKKKKKKKKEKKKKKKKKKK KKK 


SUBROUTINE eigen(c,d,dz,nz,efun,m,k,kc,w,z,a, 


& 
& 


depth, work,mstor,ifreq,isize, 
itrap, VERBOSE) 


logical ON, OFF 

parameter (ON=1, OFF=0) 

implicit real*8 (a-h,o-z) 

real*8 k(m),k_trap,k_last,kb2,pi,depth 

integer i1,j,1k,jlast,im,itrap,all 

complex*16 kc(nz) ,z(nz,nz) 

dimension c(isize) ,d(isize) ,a(nz,nz) , work (2* (2+isize) *isize) 
dimension efun(isize,m) 
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dimension mstor (isize) 
logical VERBOSE, FLAG 
data ijob/2/ 


INITIALIZE ARRAY A 


DO 10 i=1,nz 
DO 20 j=1,nz 
a(i, 3) = sorde 
20 CONTINUE 
10 CONTINUE 


INITIALIZE WORKING ARRAY 


DO i=1, ((2+1size) *isize*2) 
work(i) = 0.0 
ENDDO 


DEFINE COMPUTATIONAL DOMAIN 


nzpl =nz +1 


PRINT HEADER 
if (ifreq.eq.1) then 


write (9,*)‘'NMODE PROGRAM INITIATED’ 

write (9,*)’TOTAL MODES EVALUATED m = ‘’,m 

write (9,*)’ dz(m)=’,dz,’ nz=',nz 

write (9, *) 

write (9,*)’sample input sound speed profile at ir=1,itr=1’ 
write (9,"(4(14,1x,£8.3))") (1z7e(az)paz=lonzpl 

write (9,*)’sample input density profile at ir=1,itr=1’ 
write (9,’ (4(14,1x,£8.3))’) (iz,d(iz),iz=1,nzpl1,2) 

endif 


DEFINE Pl 


gqQaqaaQaaan 


pi = 4.d0 * datan(1.d0) 
INPUT geometrics parameters 
frequency is in hertz*2*pi=radians/sec 
dz given inm 
dz2inv = 1.d0/dz**2 
dz2inv2 = 2.d0*dz2inv 
dz2 = dz*2.d0 
write (6,*) ‘CALCULATING A MATRIX’ 
do 11 iz=2,nzpl 
set up the operator matrix a in band symmetric mode 


ai (iz.1t.nzp1) then 
dterm = (dlog(d(1z+1)) -dlog(d(iz-1)))/dz2 
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else 
dterm = 0.d0 
endif 
ieee lteazoi)) aliz-1,4z) = dz2inv-dterm 
tee (izeger2) a(lz-1,12-2) = dz2inv+edterm 
ial a(iz-1,iz-1) = (w/ce(iz))**2 - dz2inv2 


c soft upper b.c. incorporated above 
c hard lower b.c. 


a(nz,nz-1) = dz2inv2 
Cc 
© FOR FIRST FREQUENCY OR IF VERBOSE ON, OUTPUT EIGEN MATRIX 
c 
if (ifreq.eq.1 .OR. VERBOSE) then 
write (9,*)’\n check input array A(1,1 ........ Vs ag 
write (9,*)’ 2 a 
write (9,*)’ : 
ae Se : aut 
write (9,*)’ (ie rae 5,5) \n’ 
do i=1,5 
write (9,*) (a (iipeei=1, 5) 
enddo 
write (9,*)’\n check input array A(nz-5,nz-5...... Meo), 1a) 
write (9,*)’ : ad 
write (9,*)’ : at 
write(9,*)' : - 
write (9,*) ’ (1272-5 eee nz,nz) \n’ 
do i=(nz-5) ,nz 
write (9, *) (ai, ae, Je (nz -5) nz) 
enddo 
endif 
e 
c find the k*k’s and the z’s 
e 
¢ ijob = 1, calculate k only 
c¢ ijob = 2, calculate both k and z 
c ijob = 3, calculate k, z, and performance index 


write(6,*)’SOLVING FOR EIGEN VALUES & FUNCTIONS’ 
CALL eigrf (a,nz,nz,ijob,kc,z,nz,work,ier) 


write(9,*)’\n ifreq = ',ifregq,’ w= ',w 
write(9,*) ‘\n p index =’ ,work(1) 


if (work (1) .ge.50.) then 
write (6,*)’\n\nWARNING WARNING WARNING WARNING’ 
write (6,*)’P INDEX’ ,work(1),' IS TOO LARGE’ 
write (9,*)’\n\nWARNING WARNING WARNING WARNING’ 
write (9,*)’P INDEX’ ,work(1),’ IS TOO LARGE’ 
endif 


YOU DON’T REALLY WANT TO LOOK AT THIS, IT’S BORING 


if (VERBOSE) then 


write (9,*)’complex first eigenvector’ 
do 1z2=19mz , 2 


write (9, ae ez A 


QAAAAAIAN 


Ls. 


enddo 
endif 


TEST K’S FOR TRAPPED MODES 


ATG Q-On AO 


DTERMINE K AT THE BOTTOM 
write (6,*) ’DETERMINING TRAPPED MODES’ 


cb = 1677.0 

kb2 = (w/cb) **2 

itrap = 0 

write (9,*)’\n MODE FREQ:’,w/(2*pi),’\n’ 


e DEFINE DOMAIN OF EIGEN SOLUTIONS 


if(nz.lt.m) then 
all=nz 

else 
all=m 

endif 


Cc TEST FOR TRAPPED MODES 


do 14 im=1,all 
IF (abs (dimag(kc(im))).GT.0.d0 .OR. 
& dreal (kc(im)) .LT.kb2) then 
IF (VERBOSE) THEN 
WRITE (*,*) ‘MODE at im = ’,1im,’ AT FREQ’ ,w/(2*pi), 


& ‘’ Hz IS NOT TRAPPED’ 
IF (abs (dimag(kc(im))).GT.0.d0) WRITE (*,*) 
& ‘K*2 IMAG =’ ,dimag (kc (im) ) 
IF (dreal (kec(im)).LT.0.d0) WRITE(*, *) 
& ’K*2 (NEG) =’ ,dreal (kc (im) ) 
IF (dreal (kc (im) ) .LT.kb2 .AND. dreal (kc (im) ) .GT.0.d0) WRITE (*, *) 
& ’K*2 REAL < KB°2 =’, dGreal(kc(im)),’ <’,kb2 
ENDIF 
ELSE 
itrap = itrap + 1 
mstor(itrap) = im 


k_trap = dreal (kc (im) ) 
write(9,*) ’Mode Evaluated at im =’,im, 


& ’ 1s TRAPPED’ 
write (9,*)’eigen value: Kn =’,sqrt(k_trap) 
ENDIF 
14 continue 


CII III II IOIOIGI III ICICI III III III te tek 
WRITE (9,*) ‘NUMBER OF MODES \n’,itrap 


Cc kkekkkkkkkkkkkkkkkkkhkkkkkkkek kkk keke eke kk 


if (ipram.eq.1) goto 330 


NORMALIZE EIGENFUNCTION (normalize the z’s, store in efun) 


QQ 
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im=1 
mreal=0 
jJlast=0 


do 13 WHILE (im.LE.itrap) 
k_trap=0 


IF (im.EQ.1) THEN 
do ik=1,itrap 
if (dreal (kc(mstor(ik))).GT.k_trap) then 
k_trap=dreal (kc (mstor (ik) ) ) 
j=ik 
endif 
enddo 
ELSE 
do ik=1,itrap 
if (dreal (kc (mstor(ik))) .GT.k_trap.AND. 
& dreal (kc (mstor(ik))).LT.k_last) then 
k_trap=dreal (kc (mstor (ik) ) ) 
Jeik 
endif 
enddo 
ENDIF 


COUNT REALIGNED MODES 
if (j.NE. (jlast+1)) mreal = mreal+1 
iam) = k trap 
ketast = k trap 
= 0 


orthnorm .0 


do i=1,nz 
orthnorm =orthnorm + 


& (1.0/d(i+1) )* (dreal(z(i,mstor(j)))**2 + 
& dimag(z(i,mstor(j))) **2) 
enddo 


orthnorm =dsqrt (orthnorm*dz) 
TAKE REAL PART AND STORE IN EFUN ARRAY 


do i=1,nz 
efun(i,im) = dreal(z(i,mstor(j)))/orthnorm 
enddo 
jlast=j 
13 im = im+1 


write (9,999) mreal 
g99 format (14,’ modal realignments made.’ ) 


eigenfunction alignment (refr.radial dependency) 
method:test for significant value (350) 

test for sign change in gradient 

save sign of first extremum 


set up standard of comparison 


oy 


mreal=0 
FLAG=ON 


DO 250 im=1,itrap 
ext = 7002 
iz =nz 


do 251 while (dabs (efun(iz,im)) .LT.ext.AND.iz.GT.1) 
251 iz = iz-l 


e TEST FOR ZERO ENERGY AT BOTTOM (FAULTY EIGENFUNCTIONS) 


if (dfloat (efun (ifix(depth/float (dz)),im)) .EQ.0.) then 
if (FLAG) then 
write (*,*)’\nWARNING! WARNING! \n’ 
write (9,*)’\nWARNING! WARNING! \n’ 
FLAG=OFF 
endif 
write (*,1002) im 
write (99,1002) im 
endif 


Cc ALIGN ALL EIGENFUNCTIONS TO BE NEGATIVE AT THE BOTTOM 


if (dsign(1.d0,dreal (efun(iz,im))).GT.0) then 
mreal=mreal+1 
do j=1,nz 
efun(j,im) = -efun(j,im) 
enddo 
endif 


250 CONTINUE 
write (9,1001) mreal 
seven format (i4,’ modal sign changes made.’) 
1002 format (’\tMode’,13,’ May Have a Faulty Eigen Function’) 


OUTPUT SAMPLE OF EIGEN FUNCTIONS TO MODE.SYS 
IF VERBOSE SET TO ON 


AAARAN 


if (VERBOSE) then 

write (9,1003)4,k(4), (efun(iz,4) ,iz-1,nz) 
write (9,1003)3,k(3) , (efun(iz,3) ,iz2,nz) 
BP ike (9, 1003) Gua Een Cle ee 

endif 


1003 format (’ mode’,i,’ eval=’,d11.5/' efun:’ 
e S (2x,dlil 5) /22NGxsitex dane ee) 


™~ 


330 write(6,*) ‘’\n FINISHED CALCULATING EIGENMODES’ 
return 
end 


cat ttn KKKKKKKKKKKKKEKKKKEKKKKKKKKKKKKK KKK KKK KKK KKK KK IK 


c Subroutine READP 


Cc PROFILE READER. INTERPOLATES TO FILL IN PROFILES VALUES BETWEEN 
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QAQAN 


QAAQAAN 


MQM © 


Qa) 


PNEUT'S . 


- VERTICAL STEP SIZE 

- NUMBER OF POINTS IN VERTICAL ARRAY 

NUMBER OF POINTS IN WORKING ARRAY 

MUST BE > OR = ZMAX(PROFILE) /DZ + 1.5 

- FINAL VALUES OF INTERPOLATED PROFILE ARRAY 

- WORKING ARRAY OF VALUES 

NOTE: ACTUAL PROFILE VALUES MUST BE > -100.0 
OR INITALIZATION VALUES MUST BE CHANGED 


ie aS 


Mewercececec ee CCCCeCcececececccCcCCCCeCCCCECCCCCCCECCCCCCCCCCCCCCCCCCC 


SUBROUTINE READP (DZ,M,N,AF,A, IN) 
REAL*8 A,DZ,AI,Z1I,AF,IN 

INTEGER k 

DIMENSION A(N) , AF (M) , IN(N, 2) 
minor = —1 


initialize output array 


DO 7 I=1,M 
AF(I) = 0.0 
7 CONTINUE 


INITIALIZE WORKING ARRAY 
PROFILE VALUE MUST BE LESS THAN INITALIZATION VALUE 


DO 1 I=1,N 
A(I)=-100.0 
1 CONTINUE 


k=1 
ZI=IN(k,1) 
AI=IN (k, 2) 
A(1)=AI 

I=1 .54+Z1/DZ 
A(I) =AI 


2 IF(ILAST .LT. 0) THEN 
k = k+l 
ZI=IN(k,1) 
AI=IN(k, 2) 

ENDIF 


IF(ZI .GE. (M-1)*DZ) ILAST=1.5+Z1I/Dz 


EXIT LOOP IF DONE READING PROFILE 
PROFILE MORE THAN ONE PAST THE BOTTOM 
OR HAS REACHED THE LAST POINT INPUTED 


TR( Zn oRe0se) GO TO 3 
T=1.5+Z1/DZ 


ENSURE THAT N HAS BEEN DEFINED LARGE ENOUGH 


IF (feeete N) THEN 

WRITE (*,*) ‘N DEFINED TOO SMALL: N MUST BE > NZ+2’, 
& ‘ IN INPUT PROFILE, N = ',N,’ I = ‘',I1 
ELSE 


ENDIF 
IF (ILAST.GE.0) GO TO 3 


nS 


QAAAANAANA 


ene (se QaAAN oe eae COG) e) 


QAAANN 


GO TO 2 


EDED 


SET A(M+2) EQUAL TO THE LAST VALUE IN THE PROFILE 
VALUES AT DEPTHS BELOW THE LAST DEFINED POINT 
ARE GIVEN THE SAME VALUE AS THE LAST DEFINED POINT 
UNLESS A POINT DEEPER THAN THE BOTTOM IS GIVEN 
3 IF(ILAST .LT. 0) A(M+2) =A(I) 
I=l 
J=1 
4 IT=I+1 
FIND NEXT DEFINED PROFILE VALUE (FIRST VALUE ALREADY GIVEN) 
NOTE: PROFILE VALUES MUST BE LESS THAN INITIALIZATION VALUE 
IF(A(I) .LE.-100.0)GO TO 4 
IF TWO CONSECUTIVE PROFILE POINTS ARE GIVEN, NO INTERPOLATION IS 
IF(I-J.EQ.1)GO TO 6 
INTERPOLATE IN BETWEEN THE TWO DEFINED PROFILE VALUES TO FIT NEW GRID 
DO 5 K=J+1,I-1 
A (K) =A(J) +FLOAT (K-J) * (A(I) -A(J) ) /FLOAT (I-J) 
5 CONTINUE 
RESET J VALUE TO THE LAST DEFINED POINT INTERPOLATED TO 
6 J=I 
IF (ILAST.GE.0) THEN 
TF (Jd.LT. LLAST) GO TO 4 
ELSE 
IF (gir. Mt2) Gore 4 
ENDIF 
WRITE COMPUTED ARRAY INTO OUTPUT ARRAY 
OUTPUT ARRAY MUST HAVE NO MORE THAN M VALUES: 
MUST BE SMALLER THAN WORKING ARRAY 
DO K=1,M 
AF (K) =A(K) 
ENDDO 
RETURN 
END 
MODAL EXTRACTION PROGRAM 


This program, written in Fortran F77, takes the binary 


output of the modal decomposition program, and puts it into a 


usable format for modal beamforming. The output from this 


program is in two files. The first is anmeenSCll Beure 


100 


containing the eigen values with corresponding frequencies. 


The second contains the eigen function, or mode shapes for 


each frequency stored in a direct access binary format. 
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PROGRAM: "CWextract.£" 
CREATED BY GLENN A. OMANS II 
WRITTEN 30 JUN 1992 
LAST UPDATE 18 JULY 1992 
PURPOSE: TO EXTRACT MODE DATA FOR USE IN BB BEAMFORMING 


INPUT: DIRECT ACCESS BINARY DATA CREATED BY ALLFREQ.F 
FILENAME IS USUALLY modes.dabin 


OUTPUT: TWO FILES CONTAINING THE EIGEN FUNCTION (Z) FOR 
ALL FREQUENCIES AND THE EIGEN VALUES Kn’s 


EIGEN FUNCTION OUTPUT IS A SEQUENTIAL BINARY FILE 
EIGEN VALUE OUTPUT IS AN ASCII FILE 


+ + + + + + + F F FF F F OF F OF F OF OF 
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* LISTING OF VARIABLES, ARRAYS & PARAMETERS 

* 

a Dut f (array) vector holding input eigen functions for a specific 
* frequency 

2 (cha (var) frequency increment 

Zodz (var) depth increment 

* efun (array) depth by freq. matrix holding eigen function vectors 
* eval (array) elgen value vector for all frequencies & some mode # 
* EX (log) switch defining a files existance 

* feo (var) carrier/center frequency, Not used 

* fmax (var) maximum frequency in band 

* fmin (var) minimum frequency in band 

* fnameil (char) Holds binary input filename 

* fname2 (char) Holds eigen function filename 

* fname3 (char) Holds eigen value filename 

* ifreg (var) frequency counter 

* irec (var) Points to record being read from 

* itrap (var) Total number of trapped modes for a single frequency 
* izmax (var) Maximum Number of Depth Points; the min(zmax/dz and 
nzpl1-1) 

* jump (var) from input file, defines number of modes stored 

* kn (var) SORT (Eigenvalue) for specific mode & frequency 

* m (par) Maximum Number of eigen values ie. modes 

* [Must be larger then stored modes] 

* mi (var) Total number of solutions evaluated by eigen solver 
* mintrap (par) Defines Maximinimum for trapped/saved modes 

* mode (var) Input from STD IO defining mode # to extract 

* on (par) Maximum Eigen Function Length 

* [Must be larger then stored frequencies] 

* nf (var) Total number of frequencies in file 

* Nic (var) Defines the "ifreq" for the Center Frequency 

~ pi (var) PI itself in double precision 

* rol (var) Record Length, defined by input file 


al (G yal. 


* % % 


2a SQ Oe 


QAAMD 


(eee. 


c 


Q 


(var) Radial Frequency, part of eigen value output 
(array) Depth at each frequency 


zmax (par) depth cutoff for eigen function 


implicit real*8 (a-h,o-z) 
parameter (zmax=420) 

parameter (n=210,m=50) 
character*20 fnamel, fname2, fname3 


real*8 w,fc,dz,df,kn,efun, z 

real*8 pil 

real*4 fmin, fmax 

integer jump,rcl,irec,izmax 

LOGICAL EX 

DIMENSION eval (m) ,efun(n+1,m) , buff (n) 


DECLARE PI IN DOUBLE PRECISION 
pi = 4.d0 * datan(1.d0) 
PRINT HEADER 


write (*,*)’\n\n\n’ 

write (*,*)’\n\t\t\t\t\t\tPROGRAM WRITTEN: 3 AUGUST 1992’ 
write (*;*)”" \E\ENE\E Vee BY: Glenn A. Omans II’ 
write (*,*)’\n\n’ 

write (*,*)’\t*** EXTRACT EIGEN FUNCTIONS AND ’, 

& ‘VALUES FOR BB BEAMFORMING ***! 

write (*,*)’\n\nNOTE:’ 

write (*,*)’INPUT FILE IS DIRECT ACCESS BINARY’, 

& ‘\n\tUSUALLY CALLED: modes.dabin’ 

write (*,*)’\t\tCREATED FROM MODAL DECOMPOSTION PROGRAM’ 
write (*,*)’\t\t\tCALLED allfreq.f’ 

write (*,*)’\n\n\n’ 


READ MODE NUMBER FROM OPERATOR 


write(*,*)’\n INPUT MODE NUMBER TO BE EXTRACTED \n’ 
read(*,*) mode 


READ IN FILE NAMES FROM KEYBOARD 


write (*,*)’\nENTER NAME OF BINARY INPUT FILE: \n’ 
read(*,*) fnamel 

INQUIRE (FILE=fnamel , EXIST=EXx) 

IF(.NOT.EX) GOTO 2001 

write (*,*)’\nENTER NAME OF EIGEN FUNCTION OUTPUT FILE: \n’ 
read(*,*) fname2 

write (*,*)’\nENTER NAME OF EIGEN VALUE OUTPUT FILE: \n’ 
read(*,*) fname3 


OPEN INPUT AND OUTPUT FILES WITH ERROR MSG 


OPEN (UNIT=20, FILE=fnamel, STATUS=’ OLD’ , FORM=’UNFORMATED’ , 
is ACCESS=’ DIRECT’ , RECL=1700, ERR=2000) 
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AA 


QAAN 


Aaa 


*** ABOVE FILE WILL BE CLOSED & REOPENED ONCE CORRECT *** 
aoe RECORD LENGTH HAS BEEN OBTAINED a 


INQUIRE (FILE=fname2 , EXIST=EX) 
IF(EX) GOTO 2005 


OPEN EIGEN VALUE DATA FILE AFTER RECORD SIZE DETERMINED 
OPEN (UNIT=21,STATUS='NEW’ , FORM=’ FORMATED’ , FILE=fname3, 
& ERR=2010) 


READ IN DATA FROM INPUT FILE 


read (UNIT=20,REC=1) fc,fmin, fmax,df,nf,dz,mi,nzpl,jump,rcl 


write (*,*)’RECORD SIZE:',rcl 


? 
write (*,*)'’fcoc =’,fc 
write (*,*)’fmin ='’,fmin,’ fmax =’ , fmax 
write (*,*)’df =’,df,’' nf =',nf 
write (*,*)‘'dz =',dz,' Max Mode Eval =',mi,’ Eigen Size =’,nzpl 


CLOSE DATA FILE & REOPEN w/ CORRECT RECORD LENGTH 


CLOSE (20) 
OPEN (UNIT=20, FILE=fnamel, STATUS='OLD’ , FORM=’ UNFORMATED’ , 
& ACCESS=' DIRECT’ , RECL=rcl , ERR=2000) 


ENSURE SELECTED MODE NUMBER IS IN BOUNDS 
if (jump .LT.mi) then 
mintrap=j ump 
else mintrap=mi 
endif 
ERROR OUT ON NOT ENOUGH MODES STORED 
if (mode.GT.mintrap) GOTO 2020 
ERROR OUT IF VARIABLES ARE DIMENSIONED TOO SMALL 
if (m.LT.nf) then 
GOTO 2030 
endif 
Ve (meu. (izpl-1)) then 


GOTO 2040 
endif 


EXTRACT EIGEN VALUES 


do 100 ifreq=1 rit 


irec=(jump+2) * (ifreq-1) + 4 
read (UNIT=20,REC=irec) w,itrap, eval 


ERROR OUT ON "MODE NOT STORED" 


itrap=itrapt+l 


03 


if (mintrap.GT.itrap) mintrap=itrap 
if (mode.GT.itrap) goto 2020 


c OUTPUT EIGEN VALUE TO FILE 


kn = sqrt (eval (mode) ) 
write (21,*) w,kn 


100 CONTINUE 


LOAD EIGEN FUNCTION MATRIX 


im Cy 


do 300 1freg=i1 nf 
ilrec = (jump+2) * (1freq-1) +4+mode 
read (UNIT=20,REC=irec) buff 


do 200 iz=1, (nzpl1-1) 
efun(iz,ifreq) = buff (iz) 
200 continue 
300 continue 


ce 
C WRITE DEPTH VECTOR TO efun (iz, 0) 
e 


i=0 
Z=0 
do 350 iz=1, (nzpi-1) 
if (z.LE.zmax) then 
1=1+1 
efun (1z,0) =z 
endif 
350 z = Z+dz 


izmax=1 


OPEN EIGEN FUNCTION DATA FILE 


AAA 


OPEN (UNIT=22, FILE=fname2, ACCESS=’ DIRECT’ , RECL=4* (nf+1), 
& FORM=’ UNFORMATED’ , STATUS=’ NEW’ , ERR=2005) 


WRITE EIGEN FUNCTION MATRIX TO BINARY FILE 


QA 


CALL efunOUT (efun,temp,n,m,nf,izmax) 


GOTO 2050 


C 
& ERROR MESSAGES 


C 


2000 write(*,*)’\n ERROR OPENING INPUT FILE: ’,fnamel 
write (*,*)’\n INPUT FILE MUST BE BINARY’ 
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2001 


2005 


2010 


2020 


2030 


2040 


2045 


2050 


2060 


write (*,*)’\n CHECK RECORD SIZE’ 

GOTO 2050 

write (*,*)’\n ERROR INPUT FILE ’,fnamel,’ DOES NOT EXIST! !’ 
GOTO 2050 

write (*,*)’\n ERROR OPENING EIGEN VALUE DATA FILE: ‘',fname2 
write (*,*)’\n FILE MAY ALREADY EXIST’ 

GOTO 2050 

write (*,*)’\n ERROR OPENING EIGEN FUNCTION DATA FILE: ', fname3 
write (*,*)’\n FILE MAY ALREADY EXIST’ 

GOTO 2050 

write (*,*)’\n ERROR ERROR ERROR ERROR ERROR\n’ 

write (*,*) NUMBER OF MODES SELECTED’ ,mode, 


& ' EXCEEDS STORED MODES’ ,mintrap 


GOTO 2045 

write (*,*)’\n ERROR ERROR ERROR ERROR ERROR\n’ 
write (*,*) ’DIMENSION ERROR, m =’,m 

write (*,*)’MUST BE GREATER THEN’ , jump 

GOTO 2045 

write (*,*)’\n ERROR ERROR ERROR ERROR ERROR\n’ 
write (*,*)’DIMENSION ERROR, n =’,n 

write (*,*)’MUST BE GREATER THEN’ ,nf 

close (21, STATUS=’ DELETE’ ) 

close (22, STATUS=’ DELETE’ ) 

GOTO 2060 


close (21) 
close (22) 
close (20) 


STOP 
END 


kaekkkkkkkkkkkkkkkkkekkkkkkKeKKKKKKKKKKKKKKKKKKKKKRKKKRKKKKRKKKKRKKRKKKKKRKKEKEEK 
kkekkkkkkkkkkkkkkkkkkkkekkkekkekkkkkkkekkekkkkkkkkekkekkkkkkkkekkkkkekkkkkekkkkkk 


* SUBROUTINE efunOUT (efun, temp,n,m,nf,izmax) 


* 


* WRITE TO DIRECT ACCESS FILE 


* 


Kkkkk&kkkkkkekkekKekKkKkkKkKkKkKK KK KKK KKKRKKKKRKKKKRKRKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK 


400 


500 


SUBROUTINE efunOUT (efun, temp,n,m,nf,izmax) 


real*8 efun 

real*4 temp 

integer iz,n,m,nf,izmax 
dimension efun(n+1,m) ,temp (nf+1) 


INITIALIZE TEMP BUFFER 


de ifreq=! ,nf 
temp (iz) =0.0 
enddo 


do 500 iz=1,izmax 
do 400 ifreq=-0,nf 
temp (ifreq+1) = float (efun(iz,ifreq) ) 
continue 
write (22,rec=1z) temp 
continue 


RETURN 
END 


nes 


Ds. 


395. 


MODAL DECOMPOSITION SAMPLE INPUT FILE 


Odo 


2.0d0 


224. 
0 


do 


.0000000e+00 


8.0000000e+00 


RFP RRrRRPRPRPRFrRrPWUW WOU WW OWOONDT NANA AIAAARARARAUWUMUUU AH HH Hh Db WWW WWDD DDN HP PP 


.0000000e+01 
.2000000e+01 
.4000000e+01 
-6000000e+01 
-8000000e+01 
.0000000e+01 
.2000000e+01 
-4000000e+01 
.6000000e+01 
-8000000e+01 
.0000000e+01 
.2000000e+01 
.4000000e+01 
-6000000e+01 
-.8000000e+01 
.0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
-8000000e+01 
-0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
.8000000e+01 
.0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
.8000000e+01 
-0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
.8000000e+01 
.0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
.8000000e+01 
.0000000e+01 
.2000000e+01 
.4000000e+01 
.6000000e+01 
.8000000e+01 
.0000000e+02 
.0200000e+02 
.0400000e+02 
.0600000e+02 
.0800000e+02 
.1000000e+02 
.1200000e+02 
.1400000e+02 


PREP PHP PPP PP PHP PPP EPP PP PRP PP PPP P PRP PRP BPP PPP PPP PPP PPP PPP PPP PP PB 


.4790844e+03 
.4798710e+03 
.4798281e+03 
.4797952e+03 
-4797784e+03 
-4797665e+03 
.4798083e+03 
-4798116e+03 
-.4796737e+03 
.4787948e+03 
.4752782e+03 
-4738331e+03 
.4731201e+03 
.4714027e+03 
.4711126e+03 
.4703097e+03 
.4698482e+03 
.4690343e+03 
.4687429e+03 
-4687005e+03 
.4684512e+03 
.4682411e+03 
.4681399e+03 
-4679103e+03 
-4678015e+03 
-4676039e+03 
.4673772e+03 
.4672268e+03 
.4671120e+03 
.4669098e+03 
.4668384e+03 
.4666859e+03 
.4665955e+03 
.4665304e+03 
.4665222e+03 
.4664575e+03 
-4663630e+03 
.4662527e+03 
.4661536e+03 
.4661243e+03 
.4661243e+03 
.4661387e+03 
.4661316e+03 
.4661271e+03 
.4657917e+03 
.4657545e+03 
.4656872e+03 
-4656550e+03 
.4656384e+03 
.4655749e+03 
.4655615e+03 
.4655545e+03 
.4655444e+03 
.4654871e4+03 
.4654259e+03 


depth (m) 
dz (m) 
Center Freq (Hz) 
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NNNNNNNNNNNNNNNNNNNPRPRPRP PRP RP RP RPP RPP RPP RP PRP PP PPP RP RP RP PPP PPP PPP PP Pp PP PPP 


.1600000e+02 
.1800000e+02 
.2000000e+02 
.2200000e+02 
.2400000e+02 
.2600000e+02 
.2800000e+02 
.3000000e+02 
-3200000e+02 
.3400000e+02 
.3600000e+02 
-.3800000e+02 
-4000000e+02 
.4200000e+02 
-4400000e+02 
.4600000e+02 
.4800000e+02 
.5000000e+02 
.5200000e+02 
-.5400000e+02 
.5600000e+02 
.5800000e+02 
.6000000e+02 
.6200000e+02 
.6400000e+02 
.6600000e+02 
.6800000e+02 
.7000000e+02 
.7200000e+02 
.7400000e+02 
. 7600000e+02 
.7800000e+02 
.8000000e+02 
.8200000e+02 
.8400000e+02 
.8600000e+02 
-8800000e+02 
-9000000e+02 
.9200000e+02 
.9400000e+02 
.9600000e+02 
.9800000e+02 
.0000000e+02 
.0200000e+02 
.0400000e+02 
.0600000e+02 
.0800000e+02 
.1000000e+02 
.1200000e+02 
.1400000e+02 
.1600000e+02 
.1800000e+02 
.2000000e+02 
,2200000e+02 
.2400000e+02 
.2600000e+02 
.2800000e+02 
.3000000e+02 
.3200000e+02 
.3400000e+02 
.3600000e+02 


PRP PPP RPP PPP PEP PPP RPP PPP BPP PPP RPP PRP RP PPP BP PPP PPP RP PPP PP PPP PPP PPP PP Pp PP RE 


.4652586e+03 
.4651952e+03 
.4651795e+03 
.4650894e+03 
-.4650565e+03 
-4649848e+03 
.4648727e+03 
-4648223e+03 
-.4647187e+03 
.4646833e+03 
-4646635e+03 
-4645808e+03 
-4645149e+03 
.4644538e+03 
.4642631e+03 
.4642143e+03 
.4641319e+03 
-4640558e+03 
.4639944e+03 
-4638786e+03 
-4638155e+03 
.4637510e+03 
-4636782e+03 
.4636094e+03 
.4635256e+03 
.4633304e+03 
.4632558e+03 
.4632098e+03 
-4631657e+03 
.4629990e+03 
.4629380e+03 
-4629091e+03 
.4628036e+03 
.4627138e+03 
.4626069e+03 
-4625955e+03 
.4625605e+03 
.4624790e+03 
.4624430e+03 
.4624121e+03 
.4623541e+03 
.4622126e+03 
.4621396e+03 
.4620576e+03 
.4619457e+03 
.4618913e+03 
-4618832e+03 
.4618745e+03 
.4617301e+03 
.4616586e+03 
.4615538e+03 
.4615089e+03 
.4614916e+03 
.4614802e+03 
.4614468e+03 
.4613521e+03 
.4613124e+03 
-.4611986e+03 
.4611175e+03 
-4609844e+03 
-4607806e+03 
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.3800000e+02 
-4000000e+02 
.4200000e+02 
-4400000e+02 
-4600000e+02 
-4800000e+02 
.5000000e+02 
.5200000e+02 
-5400000e+02 
.5600000e+02 
-5800000e+02 
.6000000e+02 
.6200000e+02 
-6400000e+02 
-6600000e+02 
.6800000e+02 
.7000000e+02 
.7200000e+02 
.7400000e+02 


LOT? 
16779. 
L677. 
£67 s 
1677. 
oI a a 


oOooo°o°o 


i 


PRP PPP 


PRP RPP RPP PPP PP PP PP HPP BP 


.4606755e+03 
.4605276e+03 
-4598254e+03 
-4594292e+03 
-4591842e+03 
-4590498e+03 
.4587578e+03 
-4587161e+03 
-4587252e+03 
-4587113e+03 
-4587202e+03 
.4587130e+03 
.4586958e+03 
-4586952e+03 
-.4586926e+03 
.4583988e+03 
.4572621e+03 
.4572621e+03 
-4572621e+03 


Zz 


(m) 


108 


rho (cm/g*3) 


HO. 


ee. 


2s 
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